《数字信号处理》用matlab实现快速傅里叶变换_matlab快速傅里叶变换代码-CSDN博客

%FFT快速计算方法

close all;clear all;clc;
x=[2,1,-1,2,3];
nx=0:4;
K=128;
dw=2*pi/K;
k=floor((-K/2+0.5):(K/2-0.5));
X=x*exp(-j*dw*nx'*k);

figure('position',[800,300,700,200]);
m=1;n=3;
subplot(m,n,1);plot(k*dw,abs(X));
title('5点序列的DTFT和FFT');xlabel('\omega');ylabel('幅度响应');
Xd=fft([2,1,-1,2,3]);
hold on;
plot([0:4]*2*pi/5,abs(Xd),'.');grid on;

Xd1=fftshift(Xd);
subplot(m,n,2);plot(k*dw,abs(X));grid on;
xlabel('\omega');ylabel('幅度响应');title('FFT移位后');

hold on;
plot([-2:2]*2*pi/5,abs(Xd1),'.');
subplot(m,n,3);
plot(k*dw,angle(X));grid on;
title('FFT移位后');xlabel('\omega');ylabel('相位响应');

2.实验结果

二、用FFT法计算序列的线性卷积

计算序列x(n)=[2,1,3,2,1,5]与h(n)=[1,2,-1,-3]的线性卷积。

1.实验代码

%用FFT计算有限长序列的线性卷积和线性相关

close all;clear all;clc;
x=[2,1,3,2,1,5,1];
h=[1,2,-1,-3];
N=length(x)+length(h)-1;    %卷积输出对应长度。
n=0:N-1;
x=[x,zeros(1,N-length(x))]; %对x补零。
h=[h,zeros(1,N-length(h))]; %对h补零。
X=fft(x);H=fft(h);          %求x，h的FFT。
Y=X.*H;y=ifft(Y);           %求两序列的FFT相乘并求IFFT。

figure('position',[800,300,700,200]);
m1=1;m2=3;
subplot(m1,m2,1);stem(n,x,'.');grid on;
xlabel('n');ylabel('x(n)');title('序列x(n)');
subplot(m1,m2,2);stem(n,h,'.');grid on;
xlabel('n');ylabel('h(n)');title('序列h(n)');
subplot(m1,m2,3);stem(n,y,'.');grid on;
xlabel('n');ylabel('y(n)');title('序列y(n)');

2.实验结果

三、FFT计算有限长序列的频谱

1.实验代码

%计算序列的频谱（DTFT）

close all;clear all;clc;
c=[9,16,32,512];
T=0.4;
for i=1:4
    L=c(i);
    D=2*pi/(L*T);
    x=[ones(1,5),zeros(1,L-9),ones(1,4)];
    k=floor(-(L-1)/2:(L-1)/2);
    X=fftshift(fft(x,L));
    m1=2;m2=2;
    subplot(m1,m2,i);plot(k*D,real(X));grid on;
    xlabel('\omega(rad)');ylabel('X(e^j^\omega)');
    Str=['N= ',num2str(L)];
    title({Str});
end

2.实验结果

四、FFT算法实现无限长序列的频谱

x(n)=0.5nu(n)，求无限长序列的频谱。若需时域加倍长截断前后,同一频率处频谱的最大相对误差不大于0.5%,试求截断长度为多少,画出其频谱。设抽样间隔为T=0.4。

1.实验代码

%计算无限长序列的频谱

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('相角');

2.实验结果

五、用FFT算法计算连续时间非周期信号的频谱

1.实验代码

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

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

2.实验结果

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

实验六：设周期信号的频谱时两个冲激，即Xa(t)=cos(10t),试用DFT法分析其频谱。

1.实验代码

%用DFT法分析周期信号的频谱

clear all;clc;
N=input('N= ');T=0.05;n=1:N;    %原始数据。
D=2*pi/(N*T);                   %计算分辨率。
xa=cos(10*n*T);                 %求x(n)的DFT，移到对称位置。
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));grid on;
xlim([-20,20]);
xlabel('\omega');ylabel('|X(j\omega)|');
title(TITLE);                   %N=100,200,800,1200分别执行。

2.实验结果