实验五(信号抽样与恢复)和实验六(FFT算法的应用)

一、实验目的

学会用MATLAB实现连续信号的采样和重建 二、实验原理 1．抽样定理

若f(t)是带限信号，带宽为m, f(t)经采样后的频谱Fs()就是将f(t)的频谱

F()在频率轴上以采样频率s为间隔进行周期延拓。因此，当sm时，不会发生频

率混叠；而当 2．信号重建

经采样后得到信号fs(t)经理想低通h(t)则可得到重建信号f(t)，即：

s<m 时将发生频率混叠。

f(t)=fs(t)*h(t)

其中：fs(t)=f(t)(tnT)=f(nT)(tnT)

sssch(t)TsSa(ct)

所以：

f(t)=fs(t)*h(t)=f(nTs)(tnTs)*TscSa(ct)  =Tscf(nT)Sa[(tnT)]

scs上式表明，连续信号可以展开成抽样函数的无穷级数。

利用MATLAB中的sinc(t)tsin(t)来表示Sa(t)，有 Sa(t)sinc()，所以可以t得到在MATLAB中信号由f(nTs)重建f(t)的表达式如下：

f(t)=Tscf(nTs)sinc[c(tnTs)] 我们选取信号f(t)=Sa(t)作为被采样信号，当采样频率s=2m时，称为临界采样。我们取理想低通的截止频率c=m。下面程序实现对信号f(t)=Sa(t)的采样及由该采样

信号恢复重建Sa(t)：

例5-1 Sa(t)的临界采样及信号重构；

wm=1; %信号带宽

wc=wm; %滤波器截止频率 Ts=pi/wm; %采样间隔

ws=2*pi/Ts; %采样角频率 n=-100:100; %时域采样电数 nTs=n*Ts %时域采样点 f=sinc(nTs/pi);

Dt=0.005;t=-15:Dt:15;

fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); %信号重构 t1=-15:0.5:15; f1=sinc(t1/pi); subplot(211); stem(t1,f1); xlabel('kTs'); ylabel('f(kTs)');

title('sa(t)=sinc(t/pi)的临界采样信号'); subplot(212); plot(t,fa) xlabel('t'); ylabel('fa(t)');

title('由sa(t)=sinc(t/pi)的临界采样信号重构sa(t)'); grid;

例5-2 Sa(t)的过采样及信号重构和绝对误差分析

程序和例4-1类似，将采样间隔改成Ts=0.7*pi/wm , 滤波器截止频率该成wc=1.1*wm ， 添加一个误差函数 wm=1;

wc=1.1*wm; Ts=0.7*pi/wm; ws=2*pi/Ts; n=-100:100; nTs=n*Ts

f=sinc(nTs/pi);

Dt=0.005;t=-15:Dt:15;

fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); error=abs(fa-sinc(t/pi)); %重构信号与原信号误差 t1=-15:0.5:15; f1=sinc(t1/pi); subplot(311); stem(t1,f1); xlabel('kTs');

ylabel('f(kTs)');

title('sa(t)=sinc(t/pi)的采样信号'); subplot(312); plot(t,fa) xlabel('t'); ylabel('fa(t)');

title('由sa(t)=sinc(t/pi)的过采样信号重构sa(t)'); grid;

subplot(313); plot(t,error); xlabel('t');

ylabel('error(t)');

title('过采样信号与原信号的误差error(t)');

例5-3 Sa(t)的欠采样及信号重构和绝对误差分析

程序和例4-2类似，将采样间隔改成Ts=1.5*pi/wm , 滤波器截止频率该成wc=wm=1

三、上机实验内容

1．验证实验原理中所述的相关程序；

2．设f(t)=0.5*(1+cost)*(u(t+pi)-u(t-pi)) ，由于不是严格的频带有限信号，但其频谱大部分集中在[0，2]之间，带宽wm可根据一定的精度要求做一些近似。试根据以下两种情况用 MATLAB实现由f(t)的抽样信号fs(t)重建f(t) 并求两者误差，分析两种情况下的结果。 (1) wm=2 , wc=1.2wm , Ts=1; (2) wm=2 , wc=2 , Ts=2.5

实验六 FFT算法的应用

一、实验目的：加深对离散信号的DFT的理解及其FFT算法的运用。 二、实验原理和例子：N点序列的DFT和IDFT变换定义式如下：

X[k]x[n]Wn0N1knN

1， x[n]NX[k]Wk0N1knN

kn 利用旋转因子WNej2nkN具有周期性，可以得到快速算法（FFT）。

在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列

的DFT正、反变换。

例1 对连续的单一频率周期信号 按采样频率

采样，截取长度N分别选

N =20和N =16，观察其DFT结果的幅度谱。

解 此时离散序列

，即k=8。用MATLAB计

算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下：

k=8;

n1=[0:1:19];

xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1)

xlabel('t/T');ylabel('x(n)'); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1)

xlabel('k');ylabel('X(k)'); n2=[0:1:15];

xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2)

xlabel('t/T');ylabel('x(n)'); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2)

xlabel('k');ylabel('X(k)');

计算结果示于图2.1，(a)和(b)分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。

实验内容：

（1） 2N点实数序列

212cos(7n)cos(19n),n0,1,2,...,2N1x(n) N2N0,其它nN=64。用一个64点的复数FFT程序，一次算出X(k)DFT[x(n)]2N，并

绘出

X(k)。

实验要求：利用MATLAB编程完成计算，绘出相应图形。并与理论计算相比

较，说明实验结果的原因。

（1） 用以下代码实现可得图6-1所示的DFT图。 >> N=64;

>> n=0:2*N-1;

>> x=cos(2*pi*7*n/N)+1/2*cos(2*pi*19*n/N); >> X=fft(x,128); >> k=n;

>> stem(k,abs(X)) >> grid

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

图6-1

理论分析如下：

由欧拉公式得：x[n]cos(22127n)cos(19n) N2N222j(N7n)1j7n1j19n1jN(N19n)eNeNe) (eN222

对p[n]ej2knN，其2N点的DFT变换为：

2N1 p[m]p[n]en0j2mn2N2N1j2n(2km)2Nen0

当2km时，p[m]1ej22N(2km)2Nj2(2km)2N =0

1e 当2km时，即p[2k]2N

由此可得x[k]当k=14,38,90,114时有值其余为0(0k2N1）

x[14]x[114]64，x[38]x[90]32 与图6-1有相同的结论。

例2 考虑x(n)cos(0.48n)cos(0.52n)

（1） 取x(n)(0n10)时，求x(n)的DFT ： X（k）；

（2） 将（1）中的x(n)以补零方式使x(n)加长到0n100，求X（k）； （3）取x(n)(0n100)，求X（k）

要求画出x(n)和X（k），并比较（1）~（3）的结果。

解：首先定义DFT和IDFT function[XK]=dft(xn,N)

n=[0:1:N-1]; k=[0:1:N-1]; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk; XK=xn*WNnk;

function[xn]=idft(Xk,N)

n=[0:1:N-1]; k=[0;1:N-1]; WN=exp9(-j*pi/N); nk=n'*k; WNnk=WN.^(-nk);

xn=(Xk*WNnk)/N;

（1）x（n）的10点DFT

n=[0:1:99];

x=cos(0.48*pi*n)+cos(0.52*pi*n); n1=[0:1:9];y1=x(1:1:10); subplot(211);

stem(n1,y1);title('sign x(n),n[0,9]'); axis([0,10,-2.5,2.5]);text(10.2,-2.5,'n') Y1=dft(y1,10);

magY1=abs(Y1(1:1:6)); k1=0:1:5;

w0=pi*2/10*k1;

subplot(212);stem(w0/pi,magY1); title('在频域dtft分析'); xlabel('单位 pi')

由于采样的频率太小，在上图中无法确定x（n）由两个频率成分组成的。

（2） 以补零方式将x（n）加长到100个样本

n2=[0:1:99];

y2=[x(1:1:10) zeros(1,90)]; subplot(211); stem(n2,y2);

title('信号x(n),0<=n<=9+90zeros'); Y2=dft(y2,100);

magY2=abs(Y2(1:1:51));

k2=0:1:50;

w2=2*pi/100*k2; subplot(212);

plot(w2/pi,magY2);

title('DTFT幅度'); xlabel('频率(单位:pi)')

实验内容：用fft函数代替例2中的dft函数后，看看运行结果如何？

例3:用FFT分析信号频率成分

一被噪声污染的信号,很难看出它所包含的频率分量,如一个由50Hz和150Hz正弦信号构成的信号,受到均值为零、均方差为0.5的高斯随机信号的于扰,数据采样率fs=500Hz.通过FFT来分析其信号频率成分,用matlab实现如下:

fs=500; %采样频率fs=500Hz. t=0:1/fs:1; %采样周期为1/fs.

f=sin(2*pi*50*t)+sin(2*pi*150*t); % 产生信号f(t) subplot(3,1,1);plot(t,f);title('原始信号'); y=f+0.5*randn(1,length(t)); %加噪

subplot(3,1,2);plot(t,y);title('受噪声污染的信号'); N=256;

Y=fft(y,N); %对加噪信号进行FFT k=0:N-1; f=fs*k/N;

subplot(3,1,3);plot(f,abs(Y));title('FFT(幅度谱)');

(由频谱图可见,在50Hz和150Hz各出现很长的谱线,表明含噪信号y中含有这二个频率的信号.在350Hz和450Hz处也出现很长的谱线,这并不是说y中也含350Hz和450Hz的信号,这是由于采样信号的频谱是以采样频率fs为间隔周期出现而造成的)

注意: 当采样频率fs>2fm=2*150=300 Hz时,满足奈奎斯特抽样定理条件,不会产生频谱混迭现象.当fs<300 Hz时则会产生频谱混迭现象. F()1 T 0s>2 m

 不发生混叠

F()1 T

s<2 m

发生混叠

 0 F()1T s=2 m 临界 0

例4 傅里叶变换的频移特性 sssmms1sss1sss若F(f(t)]=F(ω),则

F[f(t)ej0t]F(ω-0)

例:设f(t)=sin(400 t),ω0=200.

fs=1000; %采样频率fs=1000Hz. t=0:1/fs:1;

y1=sin(400*pi*t);

y2=sin(400*pi*t).*exp(j*200*pi*t); N=512;

Y1=fft(y1,N); Y2=fft(y2,N);

subplot(3,1,1);plot(t,y1); k=0:N-1; f=fs*k/N;

subplot(3,1,2);plot(f,abs(Y1)); subplot(3,1,3);plot(f,abs(Y2));

三、实验内容 1．2N点实数序列

212cos(7n)cos(19n),n0,1,2,...,2N1x(n) N2N0,其它nN=64。用一个64点的复数FFT程序，一次算出

X(k)DFT[x(n)]2N，并绘出X(k)。

2．用fft函数代替例2中的dft函数后，看看运行结果如何？ 3． 设f(t)= 5sin(2*pi*30*t)+2sin(2*pi*60*t)+0.5sin(2*pi*90*t). . 对f(t)分别以fs1=300 Hz和fs2=150 Hz进行采样,然后将二个采样信号进行快速离散傅里叶变换(FFT),观察频谱图,指出是否产生频谱混迭现象.

4．将上题中的f(t)的频谱右移100Hz.