unit1.pas

傅立叶变换程序傅立叶变换程序傅立叶变换程序

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ShellApi , Series, ExtCtrls, TeeProcs, Chart,FFTs,aaa,Complexs,TeEngine;

type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    Edit1: TEdit;
    Label1: TLabel;
    Label2: TLabel;
    Edit2: TEdit;
    Label3: TLabel;
    Label4: TLabel;
    Edit3: TEdit;
    Label5: TLabel;
    Label6: TLabel;
    Edit4: TEdit;
    Label7: TLabel;
    Chart1: TChart;
    Series1: TLineSeries;
    Chart2: TChart;
    Chart3: TChart;
    Series2: TLineSeries;
    Series3: TLineSeries;
    Label8: TLabel;
    Button1: TButton;
    Chart4: TChart;
    Chart5: TChart;
    Series4: TLineSeries;
    Series5: TLineSeries;
    GroupBox2: TGroupBox;
    Label9: TLabel;
    Edit5: TEdit;
    Button2: TButton;
    Label10: TLabel;
    Edit6: TEdit;
    Edit7: TEdit;
    Label11: TLabel;
    Label12: TLabel;
    Panel1: TPanel;
    Label13: TLabel;
    Label14: TLabel;
    Label15: TLabel;
    Label16: TLabel;
    procedure Edit1KeyPress(Sender: TObject; var Key: Char);
    procedure FormCreate(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure Label14Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    procedure init();
  end;

var
  Form1: TForm1;
  Fs,N:integer;
  T,D:single;
  input,output:array of TComplex;
  HzLine:single;//过滤频率幅值

implementation

{$R *.dfm}
(*
  频率：单位时间内完成完整圆周运动的次数,用f表示,单位是赫兹,符号是Hz.
        当单位时间取1 s时,f=n.
  转速:单位时间内质点转过的圈数,常用符号n表示,单位符号是r/s
*)
procedure TForm1.init();
var
  maxHz:integer;
begin
  //获取原始信号每秒采样点数Fs
  Fs:=strtoint(Edit1.text);
  //获取原始信号的采样率T
  T:=1/Fs;
  //获取原始信号的长度N（总点数）
  N:=strtoint(Edit2.text);
  if N > 30720 then
  begin
    Edit2.text:='30720';
    N:=30720;
  end;
  //获取频谱的分辨率D，单位是2*PI/s
  D:=1/(N*T);
  Edit3.Text:=floattostr(D);
  (*
    计算最大频谱maxHz，因为幅值谱是偶函数，相位谱是奇函数
    都关于y轴对称，所以只需要计算一半数据N/2
    频率分辨率与频率点数的乘积D*(N/2)=1/(2*T)为最大频谱值
  *)
  maxHz:=round(1/(2*T));
  //最大频谱值（从0开始，所以最大频谱必须减一）
  Edit4.Text:=inttostr(maxHz-1);
end;

procedure TForm1.Edit1KeyPress(Sender: TObject; var Key: Char);

begin
  if Key=#13 then
  begin
    try
      init;
    except
      showmessage('确认输入了正确的数值!');
    end;
  end;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
  init;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  bfile:TextFile;
  i:integer;
  s:string;
begin
  assignFile(bfile,'b.txt');
  reset(bfile);
  setlength(input,N);
  setlength(output,N);
  series1.Clear;
  for i:=0 to N-1 do
  begin
    readln(bfile,s);
    input[i].Re:=strtofloat(s);
    series1.AddY(input[i].Re);
  end;

  //进行正向FFT变换
  ForwardFFT(input,output,N);

  //由于得到的幅频图和相频图都是关于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;
end;

procedure TForm1.Button2Click(Sender: TObject);
var
  j:integer;
  count:integer;
begin
  count:=0;//用来记录幅频图中过滤后的有效点数
  try
    Hzline:=strtofloat(Edit5.text);
  except
    showmessage('输入过滤频率数值有误');
    abort;
  end;

  for j:=0 to trunc((N-1)/2) do
  begin
    if ComplexMag(output[j])<Hzline then
    begin
      output[j].Re:=0;
      output[j].Im:=0;
      inc(count);
    end;
    series4.AddXY(j*D,ComplexMag(output[j]));
  end;

  Edit6.Text:=inttostr(trunc((N-1)/2)+1);
  Edit7.Text:=inttostr(trunc((N-1)/2)+1-count);


  InverseFFT(output,input,N);
  for j:=0 to N-1 do
  begin
    series5.AddXY(j,input[j].Re);
  end;


end;

procedure TForm1.Label14Click(Sender: TObject);
begin
    ShellExecute(Handle, 'OPEN',
    PChar('http://www.zyea.com'), nil, nil, sw_shownormal);
end;

end.