傅里叶变换的应用,matlab程序,C语言程序

1 利用FFT 计算连续时间信号的傅里叶变换

设()x t 是连续时间信号，并假设0t <时()0x t =，则其傅里叶变换由下式给出

0()()i t X x t e dt ωω∞

-=⎰ 令Γ是一个固定的正实数，N 是一个固定的正整数。当,0,1,2,,1k k N ω=Γ=-L 时，利用FFT 算法可计算()X ω。

已知一个固定的时间间隔T ，选择T 足够小，使得每一个T 秒的间隔(1)nT t n T ≤<+内，()x t 的变化很小，则式中积分可近似为

(1)0

()()()n T iwt nT n X e dt x nT ω∞+-==∑⎰

(1)01[

]()i t t n T t nT n e x nT i ωω

∞-=+==-=∑ 0

1()i T

i nT n e e

x nT i ωωω-∞-=-=∑ （27） 假设N 足够大，对于所有n N ≥的整数，幅值()x nT 很小，则式（27）变为 1

01()()i T N i nT n e X e x nT i ωωωω---=-=∑ （28）

当2/k NT ωπ=时，式（28）两边的值为 2/2/12/0211()()[]2/2/i k N

i k N N i nk N n k e e X e

x nT X k NT i k NT i k NT

ππππππ----=--==∑ （29） 其中[]X k 代表抽样信号[]()x n x nT =的N 点DFT 。最后令2/NT πΓ=，则上式变为 2/1()[]0,1,2,,12/i k N

e X k X k k N i k NT

ππ--Γ==-L （30） 首先用FFT 算法求出[]X k ，然后可用上式求出0,1,2,,1k N =-L 时的()X k Γ。 应该强调的是，式（28）只是一个近似表示，计算得到的()X ω只是一个近似值。通过取更小的抽样间隔T ，或者增加点数N ，可以得到更精确的值。如果B ω>时，幅度谱()X ω很小，对应于奈奎斯特抽样频率2s B ω=，抽样间隔T 选择/B π比较合适。如果已知信号只在时间区间10t t ≤≤内存在，可以通过对1nT t >时的抽样信号[]()x n x nT =补零，使N 足够大。

例1 利用FFT 计算傅里叶变换

如图12所示的信号

102()0t t x t -≤<⎧=⎨⎩

其它 其傅里叶变换为：

2

cos()sin()()2i X e i ωωωωωω--= 利用下面的命令，可得到()x t 的近似值和准确值。 图12 连续时间信号x(t)

N=input('Input N:');

T=input('Input T:');

%计算X(w)近似值

t=0:T:2;

x=[t-1 zeros(1,N-length(t))];

X=fft(x);

gamma=2*pi/(N*T);

k=0:10/gamma;

Xapp=(1-exp(-i*k*gamma*T))/(i*k*gamma)*X;

%计算真实值X(w)

w=::10;

Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w))./(w.*w);

plot(k*gamma,abs(Xapp(1:length(k))),'o',w,abs(Xact));

legend('近似值','真实值');

xlabel('频率(rad/s)');ylabel('|X|')

运行程序后输入N=128，T=，此时0.4909Γ=，得到实际的和近似的傅里叶变换的幅度谱如图13所示，此时近似值已经相当准确。通过增加NT 可以增加更多的细节，减少T 使得到的值更精确。再次运行程序后输入N=512，T=，此时0.2454Γ=，得到实际的和近似的傅里叶变换的幅度谱如图14所示。

图13 N=128，T=时的幅度谱

图14 N=512，T=时的幅度谱

2 利用FFT 计算离散信号的线性卷积

已知两个离散时间信号[](0,1,2,1)x n n M =-L 与[](0,1,2,1)y n n N =-L ，取L =1M N +-，对[]x n 和[]y n 右端补零，使得

[]0,1,2,,1x n n M M L ==++-L

[]0,1,2,,1y n n N N L ==++-L （31） 利用FFT 算法可以求得[]x n 和[]y n 的L 点DFT ，分别是[]X k 和[]Y k ，利用DTFT 卷积性质，卷积[]*[]x n y n 等于乘积[][]X k Y k 的L 点DFT 反变换，这也可以通过FFT 算法得到。

例2 利用FFT 计算线性卷积

已知[]0.8[]n

x n u n =，其中[]u n 为单位阶跃序列，信号[]y n 如图15所示。由于当16n >时，[]x n 很小，故M 可以取为17；N 取10，126L M N =+-=。

利用下面的Matlab 命令，可得到[]x n 、[]y n 的卷积图形如图15所示。 subplot(3,1,1);

n=0:16;

x=.^n;

stem(n,x);

xlabel('n');ylabel('x[n]');

subplot(3,1,2);

n=0:15;

y=[ones(1,10) zeros(1,6)];

stem(n,y);

xlabel('n');ylabel('y[n]')

subplot(3,1,3);

L=26;n=0:L-1;

X=fft(x,L);Y=fft(y,L);

Z=X.*Y;z=ifft(Z,L);

stem(n,z); xlabel('n');ylabel('z[n]')

图15 信号x[n]、y[n]及其卷积z[n]=x[n]*y[n]

利用下面的Matlab 命令，可得到信号x[n]、y[n]的幅度谱与相位谱如图16所示。 subplot(2,2,1);

L=26;k=0:L-1;

n=0:16;x=.^n;X=fft(x,L);

stem(k,abs(X));

axis([0 25 0 5]);

xlabel('k');ylabel('|X[k]|')

subplot(2,2,2);

stem(k,angle(X));