《数字信号处理教程》用Matlab实现快速傅里叶变换

 %用FFT计算连续时间周期信号的频谱

    
  
    
 clear all;clc;
    
  
    
 for j=1:5
    
     N=input('N=');T=0.05;n=1:N;         %原始数据。
    
     D=2*pi/(N*T);                       %计算分辨率。
    
     xa=cos(10*n*T);                     %有限长余弦序列。
    
     Xa=T*fftshift(fft(xa,N));Xa(1);     %求x(n)的DFT，移到对称位置。
    
     k=floor(-(N-1)/2:(N-1)/2);          %对于w=0对称的萘奎斯特频率下标向量。
    
     TITLE=sprintf('N=%i,L=%i',N,N*T);   %变数值为格式控制下的字符串。
    
     figure;
    
     plot(k*D,abs(Xa));
    
     xlabel('\omega');ylabel('|X(j\omega)|');
    
     title(TITLE);
    
     grid on;                            %N=100,200,800,1200分别执行。
    
 end
    
 %第五次用来清除所有生成的图形，此时N随便输入一个值即可
    
 close all;

 %用FFT计算连续周期序列的DFS并作图

    
 clear all;close all;clc;
    
  
    
 x=[ones(1,5),zeros(1,5)];
    
 N1=10;N=length(x);m=3;      %取第三个周期画图。
    
 nx=0:m*N-1;
    
 xn=x(mod(nx,N)+1);          %形成周期序列。
    
 k=[0:29]+eps;
    
 Xk=exp(-j*2*pi*k/5).*sin(pi*k/2).\sin(pi*k/10);
    
 magXk=abs(Xk);
    
 angXk=angle(Xk).*180./pi;
    
  
    
 figure;
    
 m1=4;m2=1;
    
 subplot(m1,m2,1);stem(xn,'m.','markersize',10);
    
 ylabel('x(n)');
    
 subplot(m1,m2,2);stem(k,magXk,'rP','markersize',5);
    
 ylabel('|X10(k)|');
    
 subplot(m1,m2,3);stem(k,angXk,'b*','markersize',5);
    
 ylabel('|X10(k)|');
    
 subplot(m1,m2,4);stem(k,angXk,'kx','markersize',5);
    
 xlabel('k');ylabel('argX10(k)');

 %用FFT计算序列的IDTFT

    
  
    
 clear all;close all;clc;
    
  
    
 T=1;c=[4,8];                                    %假设两种N值
    
  
    
 for i=1:2                                       %做两次循环
    
     N=c(i);D=2*pi/N;                            %取N，并求模拟频率分辨率。
    
     kl=floor(-(N-1)/2:-1/2);                    %负频率下标向量。
    
     kh=floor(0:(N-1)/2);                        %正频率下标向量。
    
     w=[kh,kl]*D;                                %将负频率移到正频率的右方。
    
     X=3+exp(-j*w)-2*exp(-j*3*w)+exp(-j*4*w);    %按新的频率排序，输入数字频率。
    
     x=ifft(X);                                  %求IFFT
    
     
    
     m1=1;m2=2;
    
     subplot(m1,m2,i);stem(T*[0:N-1],x,'rp','markersize',8);
    
     grid on;
    
     xlabel('n');ylabel('x(n)');
    
     str=['N=',num2str(N)];
    
     title(str);
    
 end

 %计算无限长序列的频谱

    
  
    
 close all;clear all;clc;
    
 T=0.4;
    
 r=1;
    
 beta=5e-3;b=0.01;
    
  
    
 while b>beta
    
     N1=2^r;n1=0:N1-1;x1=0.5.^n1;X1=fft(x1);
    
     N2=2*N1;n2=0:N2-1;x2=0.5.^n2;X2=fft(x2);
    
     k1=0:N1/2-1;k2=2*k1;                    %确定两序列同一角频率的下标。
    
     d=max(abs(X1(k1+1)-X2(k2+1)));          %对应的同一频率点上的FFT的误差的最大绝对值。
    
     X1m=max(abs(X1(k1+1)));                 %X1幅度的最大值。
    
     b=d/X1m;                                %最大相对误差的百分数。
    
     r=r+1;                                  %序列长度加倍。
    
 end
    
  
    
 k=floor(-(N2-1)/2:(N2-1)/2);                %那奎斯频率范围。
    
 D=2*pi/(N2*T);                              %角频率间隔。
    
 m1=1;m2=2;
    
 subplot(m1,m2,1);plot(k*D,abs(fftshift(X2)));grid on;
    
 title('相位普');xlabel('模拟角频率(rad/s)');ylabel('相角');
    
 subplot(m1,m2,2);plot(k*D,angle(fftshift(X2)));grid on;
    
 title('相位普');xlabel('模拟角频率(rad/s)');ylabel('相角');

 %计算连续非周期信号的频谱。

    
  
    
 close all;clear all;clc;
    
 T0=[0.05,0.02,0.01,0.01];               %4种抽样间隔。
    
 L0=[10,10,10,20];                       %4种信号记录长度，N=L0(i)/T0(i)。
    
 for i=1:4
    
    T=T0(i);N=L0(i)/T0(i);               %按顺序选用T和L。
    
    D=2*pi/(N*T);                        %频率分辨率。
    
    n=0:N-1;
    
    x=exp(-0.02*n*T).*cos(6*pi*n*T)+2*cos(14*pi*n*T);%序列。
    
    k=floor(-(N-1)/2:(N-1)/2);           %奈斯特下标向量。
    
    X=T*fftshift(fft(x));                %求x的FFT并移到对称位置。
    
    [i,X(i)];                            %检查四次循环在奈斯特频率出的幅度。
    
    m1=2;m2=2;
    
    subplot(m1,m2,i);plot(k*D,abs(X));grid on;
    
    xlabel('模拟角频率(rad/s)');ylabel('幅度');
    
    str=['T=',num2str(T),'N=',num2str(N)];
    
    title(str);                          %标题显示抽样间隔及FFT点数N。
    
 end

 %用FFT计算低通滤波器的频率响应的IDTFT

    
  
    
 clear all;close all;clc;
    
 wc=5;T=0.1*pi/wc;                           %选T为抽样间隔的1/10。
    
 N=100*2*pi/wc/T;                            %输入的杨点数。
    
 D=2*pi/(N*T);w=[0:N-1]*D;                   %模拟角频率的分辨率及频率向量。
    
 M=floor(wc/D);                              %有效频段的边界下标。
    
 H=[ones(1,M+1),zeros(1,floor(N-2*M-1)),ones(1,M)]; %按新的频率排序后的频谱。
    
 h=ifft(H/T);                                %用IFFT求h(n)/T。
    
 plot([0:N-1]*T,h);                          %用plot直线样值画出单位抽样响应。
    
 xlabel('t');ylabel('h(t)');