matlab编写关于离散的傅里叶变换

1.傅里叶级数和傅里叶变换

描述周期现象最简单的周期函数是物理学上所说的谐波函数,它是由正弦或者余弦函数来表示y(t) = Acos(wt + φ)利用三角公式，上式可以写成y(t) = Acosφcoswt 一 Asinφsint，由于φ是常数,令a= Acosφ,b= 一Asinφ,那么可以得到y(t) = acoswt十bsinwt，式中A=√a²+b²，φ=arctan(-a/b)，从这里可以看出:一个带初相位的余弦函数可以看成是一个不带初相位的正弦函数与一个不带初相位的余弦函数的合成。

谐波函数是周期函数中最简单的函数,它描述的也是最简单的周期现象，在实际中所碰到的周期现象往往要比它复杂很多,但这些复杂的函数都可以近似分解成不同频率的正弦函数和余弦函数。

这里理论知识可以参考信号与系统，那里面有对傅里叶级数和傅里叶变换的详细公式和说明。

2.离散的傅里叶变换

在信号处理中,遇到的常常不是一个函数,而是一个离散的数列,举个例子，等间隔时间取样的时间序列{x0,x1…xN-1}小.在这里数据的个数是N,一般取N为偶数，如果取2的对数对应的偶数能够加快计算速度。

下面对取值范围进行改造，首先，得到的数字信号只能在正的时间段取值,在负的时间段不能取值,但由于取的是无限长的周期序列，周期为2l，因此,把取值范围(-l,l)修改为(0,2l),这样就可以避免在负的时间段取值。

例1 计算序列x(n)的DFT.

clear all;
t= linspace(1e-3,100e-3,10);
xn = sin(100* 2*pi*t);   %产生有限序列x(n)
N= length(xn);            %获得序列的长度
WNnk= dftmtx(N);
Xk= xn*WNnk;              %计算x(n)的DFT
subplot(1,2,1);stem(1:N,xn);
subplot(1,2,2);stem(1:N,abs(Xk));

例2 已知复正弦序列x1(n)=e^jpi/8R(n),余弦序列x2(n) = cos(pi/8 n)R(n)，分别对序列求当N=16和N=8时的DFT，并绘出幅频特性曲线，并分析两种N值下DFT是否有差别及差别产生的原因。

N= 16;N1= 8;
n=0:N- 1;k=0:N1 - 1;
x1n = exp(j*pi*n/8);    %产生x1(n)
X1 = fft(x1n,N);        %计算N点DFT[x1(n)]
X2 = fft(x1n, N1);      %计算N1点DFT[x1(n)]
x2n= cos(pi*n/8);       %产生x2(n)
X3= fft(x2n, N);        %计算N点DFT[x2(n)]
X4= fft(x2n, N1);       %计算N1点DFT[x2(n)]
subplot(2,2,1);
stem(n, abs(X1),'.' );
axis([0,20,0,20]);
ylabel('|X1(k)|')
title('16点的DFT[x1(n)]')
subplot(2,2,2);
stem(n,abs(X3),'.');
axis([0,20,0,20]);
ylabel('|X2(k)|')
title('16点的DFT[x2(n)]')
subplot(2,2,3);
stem(k, abs(X2),'.');
axis([0,20,0,20]);
ylabel('|X1(k)|')
title('8点的DFT[x1(n)]')
subplot(2,2,4);
stem(k, abs(X4),'.');
axis([0,20,0,20]);
ylabel('|X2(k)|')
title('8点的DFT[x2(n)]')

3.频率分辨率

增加分辨率只有通过加大取样点N,但不是补零的方式来增加N,因为补零不是原始信号的有效信号。

例3 已知模拟信号X(k),分别取采样频率为5000Hz和1000Hz时,绘出其
傅里叶变换图

Dt= 0.0005;t=-0.005;Dt:0.005;   %模拟信号
xa= exp(-2000*abs(t));
Ts=0.0002;n=-25:1:25;           %离散时间信号
x= exp(-1000 * abs(n* Ts));
K= 500;k=0:1:K;w=pi* k/K;       
%离散时间傅里叶变换
X= x*exp(-j*n'*w);X= real(X);
w=[-fliplr(w), w(2:501)];
X=[fliplr(X),X(2:501)];
figure
subplot(2,2,1);
plot(t* 1000,xa,'.');
ylabel( 'x1(t)'); xlabel('t');
title ('离散信号');
hold on
stem(n* Ts * 1000,x);hold off
subplot(2,2,2);
plot(w/pi,X,'.');
ylabel('X1(jw)'); xlabel('f');
title('离散时间傅里叶变换');
Ts= 0.001;n=-25:1:25;
%离散时间信号
x= exp(-1000 * abs(n* Ts));
K= 500;k=0:1:K;w=pi * k/K;
%离散时间傅里叶变换
X= x*exp(-j*n'*w);X= real(X);
w=[-fliplr(w), w(2:501)];
X=[fliplr(X),X(2:501)];
subplot(2,2,3);
plot(t*1000,xa,'.');
ylabel('x2(t)'); xlabel('t');
title ('离散信号');
hold on
stem(n* Ts * 1000, x);hold off
subplot(2,2,4);
plot(w/pi,X,'.' );
ylabel( 'X2(jw)'); xlabel('f');
title('离散时间傅里叶变换');

例4 对例3中产生的离散序列x1(n)和x2(n),采用sinc函数进行内插重构。

程序如下:

clear all;
Ts1 = 0.0002;Fs1 = 1/Ts1; n1 =- 25:1:25;nTs1 =n1*Ts1;   %离散时间信号
x1= exp(-1000 * abs(nTs1));
Ts2= 0.001; 
Fs2=1/Ts2;
n2= -5:1:5;
nTs2= n2*Ts2;
x2=exp(-2000*abs(nTs2));
Dt= 0.00005;t = -0.005:Dt:0.005;                       %模拟信号重构
xa1=x1*sinc(Fs1*(ones(length(nTs1),1) * t - nTs1'*ones(1,length(t))));
xa2=x2*sinc(Fs2*(ones(length(nTs2),1) * t - nTs2'*ones(1,length(t))));
subplot(2,1,1);
plot(t*1000,xa1,'.');
ylabel('x1(t)'); xlabel('t');
title('从x1(n)重构模拟信号x1(t)');
hold on
stem(n1 * Ts1*1000,x1);
hold off
subplot(2,1,2);
plot(t*1000,xa2,'. ');
ylabel('x2(t)'); xlabel('t');
title('从x2(n)重构模拟信号x2(t)');
hold on
stem(n2*Ts2*1000,x2);
hold off

如果有网友想学习，或者对这段代码中有不明白的，可以私信我，上述的matlab代码需要仔细消化。关键是要先理解原理，傅里叶离散变换的原理也比较难懂，这里只是列出实际的例子，其更深刻的原理变换是难点，然后根据案例去理解代码。