第一句子网 - 唯美句子、句子迷、好句子大全
第一句子网 > 《数字信号处理》用matlab实现快速傅里叶变换

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

时间:2018-09-14 19:18:40

相关推荐

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

目录

用DTFT的矩阵表示法计算序列的DFT;用FFT算法计算序列的线性卷积;用FFT算法计算有限(无限)长序列的频谱;用FFT计算连续时间周期(非周期)信号的频谱。

一、用DTFT的矩阵表示法计算序列的DFT

已知x(n)=[2,1,-1,2,3],用矩阵表示法求x(ejω)=DTFT[x(n)]及x(n)的5点DFT,并将DFT的(0,2pi)范围移到与DTFT的[-π,π]重叠。

1.实验代码

%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:4L=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>betaN1=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; %序列长度加倍。endk=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:4T=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.实验结果

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。