my_fft.m

用matlab编写的关于设计滤波因子的源代码

%傅立叶法
%本方法的作用是：对指定信号进行频谱分析
%
%   输入参数：singal 为待分析的信号；输入应为256个值的数据 
%            当ISI给-1时，做时间抽选FFT正变换；当ISI给1时，做频率抽选FFT正变换
%   输出参数：H 为目标函数幅值；f 为频率，
%


%%%%%%%%%%%%%%%  傅立叶变换  %%%%%%%%%%%%%%%%%
filein='D:\data_signal';
ISI = -1;
FID = fopen('D:\data_signal.txt','r')
    [ N ] = fscanf(FID,'%d\n',1);
    [f,count] = fscanf(FID,'%f  %f\n',[2,N]);
    f=f';
fclose(FID);  

% %%%%%%%%%%%%%%  反傅立叶变换  %%%%%%%%%%%%%%%%%
% filein='D:\data_spectrum';
% ISI = 1;
% FID = fopen('D:\data_spectrum.txt','r')
%     [ N ] = fscanf(FID,'%d\n',1);
%     [f,count] = fscanf(FID,'%f  %f\n',[2,N]);
%     f=f';
% fclose(FID); 

%%%%%%%%%%%%%%%  基本参数设置  %%%%%%%%%%%%%%%%%
MMAX=1;
NN=1;
J=1;

%%%%%%%%%%%%%%%%  判断数据长度是否是2的整幂数  %%%%%%%%%%%%%%%%
for I=1:15
    M=I;
    NN=NN*2;
    if (NN==N)
        break;
    end
end

%%%%%%%%%%%%%%%%  得到输入数据的实部和虚部  %%%%%%%%%%%%%%%%
  for I=1:N
      data1(I)=f(I,1);                     %数据的实部
      data2(I)=f(I,2);                     %数据的虚部
      if ISI == -1
          signal(I)=data1(I);  %输入数据的模
      end
      if ISI == 1
         signal(I)=sqrt(data1(I).^2+data2(I).^2);  %输入数据的模
      end
  end
  
%%%%%%%%%%%%%%%%  数据的重排序  %%%%%%%%%%%%%%%%
for I=1:N
    if I<=J
        temp1=data1(J);
        temp2=data2(J);
        data1(J)=data1(I);
        data2(J)=data2(I);
        data1(I)=temp1;
        data2(I)=temp2;
    end
        M=N/2;    
        while  J>M  
            J=J-M;
            M=fix((M+1)/2);           
        end
            J=J+M;       
end

%%%%%%%%%%%%%%%%  主体计算  %%%%%%%%%%%%%%%%
while MMAX<N
    ISTEP=2*MMAX;
    for M=1:MMAX
        theta=pi*(ISI*(M-1))/MMAX;
        W1=cos(theta);
        W2=sin(theta);
        for I=M:ISTEP:N
            J=I+MMAX;
            temp1=W1*data1(J)-W2*data2(J);
            temp2=W1*data2(J)+W2*data1(J);
            data1(J)=data1(I)-temp1;
            data2(J)=data2(I)-temp2;
            data1(I)=data1(I)+temp1;
            data2(I)=data2(I)+temp2;
        end
    end
    MMAX=ISTEP;
end

%%%%%%%%%%%%%%%%  反傅立叶变换时  %%%%%%%%%%%%%%%%
if ISI>0
    for I=1:N
        data1(I)=data1(I)/N;
        data2(I)=data2(I)/N;
    end    
end

%%%%%%%%%  将正变换后的数据输出为文本格式  %%%%%%%%%%%%%%
if ISI == -1
 FID = fopen('D:\data_spectrum.txt','w')
    fprintf(FID,'%20d \n',N);
    for I=1:N
        fprintf(FID,'%20.10f  %20.10f\n',data1(I),data2(I));
    end 
 fclose(FID);  
end

 %%%%%%%%%  将反变换后的数据输出为文本格式  %%%%%%%%%%%%%%
if ISI == 1
  FID = fopen('D:\data_signal.txt','w')
     fprintf(FID,'%20d \n',N);
     for I=1:N
         fprintf(FID,'%20.10f  %20.10f\n',data1(I),data2(I));
     end 
  fclose(FID); 
end 

%%%%%%%%%%%%%%%%  求出振幅和相位  %%%%%%%%%%%%%%%%
for I=1:N
    H(I)=complex(data1(I),data2(I));           %  振幅的计算 
    module(I)=sqrt(data1(I)^2+data2(I)^2);     %  模的计算
    P(I)=atan(data2(I)/data1(I));              %  相位的计算
end
%%%%%%%%%%%%%%%%%%%%%%  画图  %%%%%%%%%%%%%%%%%%%%%%
if ISI == -1
    figure(1);
    subplot(311); %图形显示分割窗口 
    plot(signal);
    title('输入信号');
    subplot(312); %图形显示分割窗口 
    plot(module);
    title('频谱分析振幅谱');
    subplot(313); %图形显示分割窗口 
    plot(P);
    title('频谱分析相位谱');
end
if ISI == 1
   figure(2);
   subplot(211); %图形显示分割窗口 
   plot(signal);
   title('输入信号');
   subplot(212); %图形显示分割窗口 
   plot(module);
   title('反傅立叶变换后振幅谱');     
end