кореляції
HARMONI, // кількість додатніх гарм. LSAMPLES/2
f1,f2; // смуга частот для анал_зу
FILE *wavfile;
RIFFHead rh;
PCMWAVEFORMAT FMTChunk;
DATAHead dh;
unsigned long int dataskip, i, j, NRead;
double sred=0, betta;
if (OpenDialog1->Execute())
{
Label5->Caption = StrRScan(OpenDialog1->FileName.c_str(),'\\');
LSAMPLES = StrToInt(ComboBox1->Text);
HARMONI = LSAMPLES>>1;
NSAMPLES = 3*(LSAMPLES>>2);
MSAMPLES = LSAMPLES>>2;
betta = StrToFloat(Edit4->Text); // коеф_ц_єнт для ф_льтра
// выделяем память под массивы
in_r = new short int[NSAMPLES];
covar = new double[LSAMPLES];
out_rel = new double[LSAMPLES];
out_im = new double[LSAMPLES];
out_mod = new double[LSAMPLES];
// проверяем наличие выделенной памяти
CHECKPOINTER(in_r);
CHECKPOINTER(covar);
CHECKPOINTER(out_rel);
CHECKPOINTER(out_im);
CHECKPOINTER(out_mod);
f1 = StrToInt(Edit2->Text);
f2 = StrToInt(Edit3->Text);
// чтение WAV-файла
if ((wavfile = fopen(OpenDialog1->FileName.c_str(),"rb")) != NULL)
{
fread(&rh,sizeof(rh),1,wavfile);
fread(&FMTChunk,sizeof(FMTChunk),1,wavfile);
fread(&dh,sizeof(dh),1,wavfile);
SamplRate = FMTChunk.wf.nSamplesPerSec;
//Label5->Caption = "SampleRate=" + IntToStr(SamplRate);
dataskip = floor(StrToFloat(Edit1->Text)*SamplRate)*2;
dataskip = (dataskip>=dh.ChunkSize) ? (dh.ChunkSize-NSAMPLES) : dataskip;
if (dataskip) fseek(wavfile,dataskip,SEEK_CUR); // пропустить сколько надо
NRead = fread(in_r,sizeof(unsigned short int),NSAMPLES,wavfile); // считать данные
if (NRead == NSAMPLES)
{
Series1->Clear();
switch (RadioGroup1->ItemIndex) {
case 0:
for(i=0;i<NSAMPLES;i++)
sred += covar[i] = in_r[i];
break;
case 1:
for(i=0;i<NSAMPLES-1;i++)
sred += covar[i] = in_r[i+1]-betta*in_r[i];
covar[i] = 0;
break;
case 2:
for(i=0;i<NSAMPLES-1;i++)
sred += covar[i] = in_r[i+1]+betta*in_r[i];
covar[i] = 0;
break;
case 3:
for(i=0;i<NSAMPLES;i++)
if (in_r[i]>0) covar[i] = 1;
else
if (in_r[i]<0) covar[i] = -1;
else covar[i] = 0;
break;
}
sred /=(float)NSAMPLES;
Label4->Caption = "Mx = " + FloatToStrF(sred,ffGeneral,7,12);
for(i=0;i<NSAMPLES;i++)
{
covar[i] = covar[i] - sred; // центруємо
Series1->AddXY((double)i/SamplRate*1000.,covar[i],"",clRed);
}
for(i=NSAMPLES;i<LSAMPLES;i++) covar[i] = 0; // решту обнуляємо
// ---------------- кореляц_я за допомогою FFT
fft_double(LSAMPLES,0,covar,NULL,out_rel,out_im);
for(i=0;i<LSAMPLES;i++)
out_mod[i] = out_rel[i]*out_rel[i]+out_im[i]*out_im[i];
fft_double(LSAMPLES,1,out_mod,NULL,out_rel,out_im);
if (CheckBox2->State) // нормування
{
for(i=1;i<MSAMPLES;i++)
out_rel[i] = out_rel[i]/out_rel[0];
out_rel[0] = 1;
}
else
for(i=0;i<MSAMPLES;i++)
out_rel[i]/=(double)(LSAMPLES);
// ----------------- згладжувальне в_кно Бартлетта
Series9->Clear();
for(i=0;i<MSAMPLES;i++)
{
if (CheckBox1->State) // згладжування даних
covar[i]=out_rel[i]*(1.-i/(double)(MSAMPLES));
else
covar[i]=out_rel[i];
Series9->AddXY((double)i/SamplRate*1000.,covar[i],"",clRed);
}
for(i=MSAMPLES;i<LSAMPLES;i++)
covar[i] = 0;
// ----------------- отримання згладженого спектру
fft_double(LSAMPLES,0,covar,NULL,out_rel,out_im);
Series17->Clear();
for(i=1;i<HARMONI;i++)
{
out_mod[i] = sqrt(out_rel[i]*out_rel[i]+out_im[i]*out_im[i]);
Series17->AddXY((double)SamplRate*i/LSAMPLES,out_mod[i],"",clRed);
}
Label13->Caption = FloatToStrF(Energy(SamplRate,HARMONI,f1,f2,out_mod)/1000000.,ffGeneral,7,12);
} // if NRead
fclose(wavfile); // закрыть файл
} // if open file
} // if opendialog
delete out_mod;
delete out_im;
delete out_rel;
delete covar;
delete in_r;
}
//--------------------------------------------------------------
void __fastcall TForm1::ButtonClear(TObject *Sender)
{
Series1->Clear(); Series9->Clear(); Series17->Clear();
Series2->Clear(); Series10->Clear(); Series18->Clear();
Series3->Clear(); Series11->Clear(); Series19->Clear();
Series4->Clear(); Series12->Clear(); Series20->Clear();
Series5->Clear(); Series13->Clear(); Series21->Clear();
Series6->Clear(); Series14->Clear(); Series22->Clear();
Series7->Clear(); Series15->Clear(); Series23->Clear();
Series8->Clear(); Series16->Clear(); Series24->Clear();
Label5->Caption = ""; Label13->Caption = "";
Label6->Caption = ""; Label14->Caption = "";
Label7->Caption = ""; Label15->Caption = "";
Label8->Caption = ""; Label16->Caption = "";
Label9->Caption = ""; Label17->Caption = "";
Label10->Caption = ""; Label18->Caption = "";
Label11->Caption = ""; Label19->Caption = "";
Label12->Caption = ""; Label20->Caption = "";
}
//--------------------------------------------------------------
ДОДАТОК Б
Текст програми Noise Flowmeter для проведення експериментальних досліджень властивостей та характеристик потоку
//--------------------------------------------------------------
#include <vcl.h>
#include <stdio.h>
#include <mmsystem.h>
#pragma hdrstop
#include "metrolog.h"
#include "fourier.h"
#include "consumption.h"
#include "consumption.cpp"
#include "FFTMISC.C"
#include "FOURIERD.C"
//--------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//--------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
: TForm(Owner)
{
}
//--------------------------------------------------------------
void __fastcall TForm1::Button1Click(TObject *Sender)
{
short int *in_r=NULL; // вказівник на вихідний масив
double *covar=NULL,
*out_rel=NULL,
*out_im=NULL,
*out_mod=NULL; // вказівник на масив гармонік
unsigned long int SamplRate; // SamplingRate = семплів / сек
unsigned int
NSAMPLES, // довжина виборки
LSAMPLES, // кількість точок для ДПФ
MSAMPLES, // довжина кореляц_ї
HARMONI, // кол-во положительных гармоник LSAMPLES/2
f1,f2; // смуга частот для аналізу
FILE *wavfile;
RIFFHead rh;
PCMWAVEFORMAT FMTChunk;
DATAHead dh;
unsigned long int dataskip, timestop, timestop1, i, j, NRead, curpos;
unsigned long int i_vidriz, i_mitteve;
double Q_vidriz, Q_vytrata, Q_tick, Q,
Qz[3],r1=0,r2=0, sred=0;
if (OpenDialog1->Execute())
{
Label5->Caption = StrRScan(OpenDialog1->FileName.c_str(),'\\');
Label6->Caption = "Середнє = ";
LSAMPLES = StrToInt(ComboBox1->Text);
HARMONI = LSAMPLES>>1;
NSAMPLES = 3*(LSAMPLES>>2);
MSAMPLES = LSAMPLES>>2;
// выделяем память под массивы
in_r = new short int[NSAMPLES];
covar = new double[LSAMPLES];
out_rel = new double[LSAMPLES];
out_im = new double[LSAMPLES];
out_mod = new double[LSAMPLES];
// проверяем наличие выделенной памяти
CHECKPOINTER(in_r);
CHECKPOINTER(covar);
CHECKPOINTER(out_rel);
CHECKPOINTER(out_im);
CHECKPOINTER(out_mod);
f1 = StrToInt(Edit3->Text);
f2 = StrToInt(Edit4->Text);
// чтение WAV-файла
if ((wavfile = fopen(OpenDialog1->FileName.c_str(),"rb")) != NULL)
{
fread(&rh,sizeof(rh),1,wavfile);
fread(&FMTChunk,sizeof(FMTChunk),1,wavfile);
fread(&dh,sizeof(dh),1,wavfile);
SamplRate = FMTChunk.wf.nSamplesPerSec;
//Label5->Caption = "SampleRate=" + IntToStr(SamplRate);
timestop = ceil(StrToFloat(Edit2->Text)*SamplRate); // пакетний
timestop1 = ceil(StrToFloat(Edit5->Text)*SamplRate); // загальний
dataskip = floor(StrToFloat(Edit1->Text)*SamplRate)*2;
dataskip = (dataskip>=dh.ChunkSize) ? (dh.ChunkSize-SamplRate*timestop) : dataskip;
timestop1 = (timestop1) ? timestop1 : (dh.ChunkSize/FMTChunk.wf.nBlockAlign);
timestop = (timestop) ? timestop : timestop1;
ProgressBar1->Position = 0;
ProgressBar1->Max = (timestop1 < dh.ChunkSize/FMTChunk.wf.nBlockAlign) ? ceil(timestop1/NSAMPLES+1) : (dh.ChunkSize/FMTChunk.wf.nBlockAlign/NSAMPLES+1);
if (dataskip) fseek(wavfile,dataskip,SEEK_CUR); // пропустить сколько надо
i_vidriz=0; Q_vytrata=0;
Q=0; j=0;
Series1->Clear();
Series2->Clear();
Series3->Clear();
do {
i_vidriz++; Q_vidriz=0; i_mitteve=0;
do {
NRead = fread(in_r,sizeof(unsigned short int),NSAMPLES,wavfile); // считать данные
if (NRead != NSAMPLES)
for (i=NRead; i<NSAMPLES; i++)
in_r[i]=0;
i_mitteve++; j++; ProgressBar1->StepIt();
switch (RadioGroup1->ItemIndex) {
case 0:
for(i=0;i<NSAMPLES;i++)
covar[i] = in_r[i];
break;
case 1:
for(i=0;i<NSAMPLES-1;i++)
covar[i] = in_r[i+1]-in_r[i];
covar[i] = 0;
break;
case 2:
for(i=0;i<NSAMPLES-1;i++)
covar[i] = in_r[i+1]+in_r[i];
covar[i] = 0;
break;
case 3:
for(i=0;i<NSAMPLES;i++)
if (in_r[i]>0) covar[i] = 1;
else
if (in_r[i]<0) covar[i] = -1;
else covar[i] = 0;
break;
}
if (CheckBox1->State) // згладжування даних
for(i=0;i<NSAMPLES;i++)
covar[i] = covar[i]*0.5*(1.-cos(2.*PI*i/NSAMPLES));
for(i=0;i<NSAMPLES;i++) sred+=covar[i]; //визначаємо середнє виб_рки
sred /=(double)NSAMPLES;
for(i=0;i<NSAMPLES;i++) covar[i] = covar[i] - sred; // центруємо
for(i=NSAMPLES;i<LSAMPLES;i++) covar[i] = 0; // решту обнуляємо
// ---------------- кореляц_я за допомогою FFT
fft_double(LSAMPLES,0,covar,NULL,out_rel,out_im);
for(i=0;i<LSAMPLES;i++)
out_mod[i] = out_rel[i]*out_rel[i]+out_im[i]*out_im[i];
fft_double(LSAMPLES,1,out_mod,NULL,out_rel,out_im);
if (CheckBox2->State) // нормування
{
for(i=1;i<MSAMPLES;i++)
out_rel[i] = out_rel[i]/out_rel[0];
out_rel[0] = 1;
}
else
for(i=0;i<MSAMPLES;i++)
out_rel[i]/=(double)(NSAMPLES);
// ----------------- згладжувальне в_кно Бартлетта
for(i=0;i<MSAMPLES;i++)
covar[i]=out_rel[i]*(1.-i/(double)(MSAMPLES))/(double)(NSAMPLES);
for(i=MSAMPLES;i<LSAMPLES;i++)
covar[i] = 0;
// ----------------- отримання згладженого спектру
fft_double(LSAMPLES,0,covar,NULL,out_rel,out_im);
for(i=0;i<LSAMPLES;i++)
out_mod[i] = sqrt(out_rel[i]*out_rel[i]+out_im[i]*out_im[i]);
// Q_tick = Energy(SamplRate,HARMONI,f1,f2,out_mod);
Q_tick = log(3.*sqrt(2.*out_mod[0]));
Q_vidriz += Q_tick;
Q += Q_tick;
Series1->AddXY((double)NSAMPLES*j/SamplRate,Q_tick,"",clRed);
} // if NRead
} while (!feof(wavfile) && j*NSAMPLES < timestop1 && i_mitteve*NSAMPLES < timestop);
Q_vytrata += Q_vidriz/i_mitteve;
// Label6->Caption = Label6->Caption + FloatToStrF(Q_vidriz/i_mitteve,ffGeneral,7,12);
// Series2->AddXY((double)NSAMPLES*(j-i_mitteve)/SamplRate,Q_vidriz/i_mitteve,"",clWhite);
// Series2->AddXY((double)NSAMPLES*j/SamplRate,Q_vidriz/i_mitteve,"",clWhite);
} while (!feof(wavfile) && j*NSAMPLES < timestop1);
Label6->Caption = "Середнє пакет_в = " + FloatToStrF(Q_vytrata/i_vidriz,ffGeneral,5,8);
Label8->Caption = "Середнє загальне = " + FloatToStrF(Q/j,ffGeneral,5,8);
fclose(wavfile); // закрыть файл
} // if open file
} // if opendialog
delete out_mod;
delete out_im;
delete out_rel;
delete covar;
delete in_r;
}
//--------------------------------------------------------------