unit1.cpp

fft变换算法的c++源代码

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "Unit1.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma link "PERFGRAP"
#pragma resource "*.dfm"
TForm1 *Form1;

int Fs,N;
float T,D;
complex <double> *input,*output;
float HzLine;

//---------------------------------------------------------------------------
void __fastcall TForm1::Init()
{
    int maxHz;
  //获取原始信号每秒采样点数Fs
  Fs=StrToInt(Edit1->Text);
  //获取原始信号的采样率T
  T=1.0/Fs;
  //获取原始信号的长度N（总点数）
  N=StrToInt(Edit2->Text);
  if(N > 30720)
  {
    Edit2->Text="30720";
    N=30720;
  }
  //获取频谱的分辨率D，单位是2*PI/s
  D=1/(N*T);
  Edit3->Text=FloatToStr(D);
/*
    计算最大频谱maxHz，因为幅值谱是偶函数，相位谱是奇函数
    都关于y轴对称，所以只需要计算一半数据N/2
    频率分辨率与频率点数的乘积D*(N/2)=1/(2*T)为最大频谱值
  */
  maxHz=(int)(1/(2*T));
  //最大频谱值（从0开始，所以最大频谱必须减一）
  Edit4->Text=IntToStr(maxHz-1);
}
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
Init();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Edit1KeyPress(TObject *Sender, char &Key)
{
 if(Key==VK_RETURN)
    try
    {
    Init();
    }
    catch(...)
    {
    ShowMessage("确认输入了正确的数值");
    }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button1Click(TObject *Sender)
{
  FILE*fp;
  int i;
  String s;
  double dd;
  TFFT FFT;
  fp=fopen("b.txt","r");
  input=new complex <double>[N];
  output=new complex <double>[N];
//  setlength(input,N);
//  setlength(output,N);
//  series1.Clear;
  Image1->Canvas->MoveTo(0,0);
  for(i=0;i<N;i++)
  {
  fscanf(fp,"%lf",&dd);
  input[i]=complex <double>(dd,0);
 Image1->Canvas->LineTo((i*Image1->Width/N),Image1->Height/2-dd*5);
 Image1->Refresh();
//  PerformanceGraph1->DataPoint(clRed,(long)(500+10*dd));
  }
//    PerformanceGraph1->Update();
//  for i:=0 to N-1 do
//  begin
//    readln(bfile,s);
//    input[i].Re:=strtofloat(s);
//    series1.AddY(input[i].Re);
//  end;

  //进行正向FFT变换
  FFT.ForwardFFT(input,output,N);
  Image2->Canvas->MoveTo(0,Image2->Height/2);
  Image3->Canvas->MoveTo(0,Image2->Height/2);
 for(i=0;i<(int)((N-1)/2);i++)
 {
 Image2->Canvas->LineTo(i*D,Image2->Height-2-0.05*sqrt(output[i].real()*output[i].real()+output[i].imag()*output[i].imag()));
 Image2->Refresh();
 Image3->Canvas->LineTo((i*D),Image2->Height/2-18*atan2(output[i].imag(),output[i].real())/pi);
 Image3->Refresh();
 }
  //由于得到的幅频图和相频图都是关于Y对称，故只需N/2的数据
  //将(N-1)/2取最大整值:trunc((N-1)/2)+1
  //然后再-1，因为起点为0
  //最后得到trunc((N-1)/2)
//  series2.Clear;
//  series3.Clear;
//  for i:=0 to trunc((N-1)/2) do
//  begin
    //频谱幅值
//    series2.AddXY(i*D,ComplexMag(output[i]));
    //频谱相位
//    series3.AddXY(i*d,ComplexPhase(output[i]));
//  end;

  Button2->Enabled=True;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button2Click(TObject *Sender)
{

  int j;
  int count;
  TFFT FFT;
  count=0;//用来记录幅频图中过滤后的有效点数
  try{
    HzLine=StrToFloat(Edit5->Text);
  }
  catch(...)
  {
    ShowMessage("输入过滤频率数值有误");
    return;
   }
  Image4->Canvas->MoveTo(0,Image4->Height/2);
  Image5->Canvas->MoveTo(0,Image5->Height/2);
  for(j=0;j<=(int)((N-1)/2);j++)
  {
   if(sqrt(output[j].real()*output[j].real()+output[j].imag()*output[j].imag())<HzLine)
    {
      output[j]=complex<double>(0,0);
      
      count++;
    }
     Image4->Canvas->LineTo(j*D,Image4->Height-2-0.05*sqrt(output[j].real()*output[j].real()+output[j].imag()*output[j].imag()));
     Image4->Refresh();
    //series4.AddXY(j*D,ComplexMag(output[j]));
  }
  Edit6->Text=IntToStr((N-1)/2+1);
  Edit7->Text=IntToStr((N-1)/2+1-count+1);


  FFT.InverseFFT(output,input,N);
  for(j=0;j<=N-1;j++)
   {
    Image5->Canvas->LineTo((j*Image5->Width/N),Image5->Height/2-5*input[j].real());
     Image5->Refresh();
   //    series5.AddXY(j,input[j].Re);
   }
}
//---------------------------------------------------------------------------