基于MATLAB的离散非周期信号频域分析

1、基于MATLAB的离散信号频域分析、快速傅里叶变换与采样定理离散信号频域分析周期离散方波信号频域分析与周期模拟信号一样，周期离散信号同样可以展开成傅里叶级数形式，并得到离散傅里叶级数（DFS）Xk=1Nn=-N2N2xne-jkn k=0,1,2,N-1 上式可以看成周期离散信号x(n)的离散傅里叶级数展开。xn=k=0N-1Xkejkn上式是DFS的反变换，记作IDFS并且称X(k)与x(n)构成一对离散傅里叶级数变换对。（以上两式中=2/N）在MTALAB中，DFS通过建立周期延拓函数语句实现：function Xk=DFS(n,x,N)if Nlength(x) n=0:N-1; x=x

2、 zeros(1,N-length(x);endk=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=x*WNnk;end建立一个离散非周期方波信号xn=RNn=1, &0nN-10, &其他R4n通过周期延拓后所得的周期序列利用DFS计算实现代码如下：clear all;close all;clc;n=0:3;x=ones(1,4);X=fft(x,1024);Xk1=DFS(n,x,4);Xk2=DFS(n,x,8);figure(1);plot(-1023:2048)/2048*8,abs(X) abs(X) abs(X),-);hold on;s

3、tem(-4:7,abs(Xk1) abs(Xk1) abs(Xk1),LineWidth,2);grid;figure(2);plot(-1023:2048)/2048*16,abs(X) abs(X) abs(X),-);hold on;stem(-8:15,abs(Xk2) abs(Xk2) abs(Xk2),LineWidth,2);grid;set(gcf,color,w);运行后得到的是分别以4和8为周期延拓后的R4n频谱：即第一幅图表示的是周期序列 xn=1 -n+ 的频谱，第二幅图表示的是周期序列xn=1, &4kn4+4k0, &4k-4n4k的频谱。两图中的包络线表示的是通

4、过快速傅里叶变换（FFT）所得到的频谱线。（二）非周期离散方波信号频域分析对于非周期离散方波信号，可采用离散时间傅里叶变换DTFT进行分析。X=n=-+x(n)e-jn上式为离散时间信号x(n)的离散时间傅里叶变换（DTFT）。xn=1202Xejnd上式为X的离散时间傅里叶反变换（IDTFT）。由于：i=-+x(i)length(x) n=0:N-1; x=x zeros(1,N-length(x);endk=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=x*WNnk;End建立一个离散非周期方波信号xn=RNn=1, &0nN-10, &其他R8

5、n的离散傅里叶变换Xej利用DFT计算实现代码如下：clear all;close all;clc;n=0:7;x=ones(1,8);X=fft(x,1024);Xk2=DFT(n,x,16);figure(1);plot(-1023:2048)/2048*32,abs(X) abs(X) abs(X),-);hold on;stem(-16:31,abs(Xk2) abs(Xk2) abs(Xk2),LineWidth,2);grid;figure(2);plot(-1023:2048)/2048*32,angle(X) angle(X) angle(X),-);hold on;stem(

6、-16:31,angle(Xk2) angle(Xk2) angle(Xk2),LineWidth,2);grid;set(gcf,color,w);运行后分别得到该离散非周期方波信号的幅频特性与相频特性：幅频特性相频特性两图中的包络线表示的是通过快速傅里叶变换（FFT）所得到的频谱线。离散傅里叶变换是傅里叶变换在时域、频域均离散化的形式，因而与其他傅里叶变换有着相似的性质。但是它又是从傅里叶级数派生而来的，所以又具有一些与其他傅里叶变换不同的特性，最主要的是圆周位移性质和圆周卷积性质。快速傅里叶变换（FFT）快速傅里叶变换，简称FFT，是计算DFT的快速算法，习惯上是指以库利和图基算法为基础

7、的一类高效算法。根据快速傅里叶变换基本思路以及基2FFT算法，在MTALAB中，FFT通过建立函数实现：function y=fft(x)m=nextpow2(x); N=2m;if length(x)0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facel

8、ighting,phong); 二、理想高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,inte

9、rp,Edgecolor,none,Facelighting,phong); Butterworth低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(

10、2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(i,j); Hd(i,j)=1/(tn+1); endendY=fft2

11、(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);D=100/si

12、ze(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolo

13、r,none,Facelighting,phong); 高斯高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(

14、Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D

15、0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形高通滤波器IA=imread(lena.bmp);IB=imrea

16、d(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2

17、,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 用其他方法编写的理想低通、理想高通、Butterworth低通、同态滤波程序理想低通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sq

18、rt(i-u0)2+(j-v0)2); if d=d0 h=1; else h=0; end y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 理想高通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshi

19、ft(k);N1,N2=size(g);n=2;d0=10;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=0; else h=1; end y(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); Butterworth低通i1=imread(lena.bmp);i2=imnoise(i1,salt