数字信号实验2

数字信号处理实验学号：110220217 姓名：齐兴
实验一离散傅里叶变换的性质
实验内容和步骤
x1[n]，x2[n] 07
n
≤≤为长度N=8的实序列，x1[n]=[1 3 5 3 6 8 3 9]，x2[n]=[2 4 3 6 7 9
0 2 ]，采用MATLAB编程验证傅里叶变换的如下性质
1.线性特性
a.给出序列x1[n]的傅里叶变换X1[k]，并画出其幅度谱和相位谱
b.给出序列x2[n]的傅里叶变换X2[k] ，并画出其幅度谱和相位谱
c.给出序列Z=2*X1[k]+6*x2[k],并与序列2*x3[n]+6*x4[n]的傅里叶变换比较，
程序清单：
clc;
clear;
x1=[1 3 5 3 6 8 3 9]; x2=[2 4 3 6 7 9 0 2]; X1=fft(x1);
X2=fft(x2);
a1=abs(X1);
a2=angle(X1);
b1=abs(X2); b2=angle(X2);
subplot(321);stem(a1);title('x1幅度谱'); subplot(322);stem(a2);title('x1相位谱'); subplot(323);stem(b1);title('x2幅度谱'); subplot(324);stem(b2);title('x2相位谱'); Z=2*X1+6*X2;
W=2*x1+6*x2;
subplot(325);stem(Z);
subplot(326);stem(W);
图像：
2.时移特性
给出序列x1[n]右移3位(循环移位)后的傅里叶变换的幅度谱和相位谱，并和原始序列的幅度谱和相位谱相比较
程序清单：
clc;
clear;
x1=[1 3 5 3 6 8 3 9]; x2=circshift(x1,[0,3]) ;
X1=fft(x1);
X2=fft(x2);
a1=abs(X1); a2=angle(X1);
b1=abs(X2);
b2=angle(X2);
subplot(221);stem(a1);title('x1原始幅度谱'); subplot(222);stem(a2);title('x1原始相位谱'); subplot(223);stem(b1);title('x1右移3位幅度谱'); subplot(224);stem(b2);title('x1右移3位相位谱');
图像：
3.对称性
（1）利用x1[n]构造圆周共轭对称序列和圆周共轭反对称序列，讨论如下问题（a）画出圆周共轭对称序列的傅里叶变换的幅度谱和相位谱
（b）画出圆周共轭对称序列的傅里叶变换的实部和虚部
程序清单：
clc;
clear;
x1=[1 3 5 3 6 8 3 9] x2=fliplr(x1)
x3=circshift(x2,[1,1]) xe=(x1+conj(x3))/2 xo=(x1-conj(x3))/2 subplot(421);stem(a1);title('圆周共轭对称序列幅度谱');
subplot(422);stem(a2);title('圆周共轭对称序列相位谱');
subplot(423);stem(a3);title('圆周共轭对称序列实部');
subplot(424);stem(a4);title('圆周共轭对称序列虚
X1=fft(xe);
a1=abs(X1);
a2=angle(X1); a3=real(X1);
a4=imag(X1);
X2=fft(xo);
b1=abs(X2);
b2=angle(X2);
b3=real(X2);
b4=imag(X2); 部'); subplot(425);stem(b1);title('圆周共轭反对称序列幅度谱'); subplot(426);stem(b2);title('圆周共轭反对称序列相位谱'); subplot(427);stem(b3);title('圆周共轭反对称序列实部'); subplot(428);stem(b4);title('圆周共轭反对称序列虚部');
图像：
（c ） 利用x1[n]构造共轭对称序列和共轭反对称序列（第二章学习），讨论与圆周共轭
对称圆周共轭反对称的区别
clc;
clear all;
close all;
x1=[1 3 5 3 6 8 3 9];
N=length(x1);
n_axis=[0:N-1];
x1N=zeros(1,N);
x1N(0+1)=x1(0+1);
for num=0:N-1
if N-num==N
index=1;
else
index=(N-num)+1;
end x1N(num+1)=x1(index);
endx1_e=1/2*(x1+x1N);
x1_o=1/2*(x1-x1N);
figure(1);
subplot(3,1,1);stem(n_axis,x1);title( ['原序列'] ); subplot(3,1,2);stem(n_axis,x1_e);title('共轭对称部分');
subplot(3,1,3);
stem(n_axis,x1_o);title('共轭反对称部分');
图像：
012345670
5
10
原序列
012345670
5
10
共轭对称部分01234567-50
5
共轭反对称部分
圆周共轭对称的主值序列是共轭对称序列，圆周共轭反对称的主值序列是共轭反对称序列；
共轭对称序列的周期延拓是圆周共轭序列，共轭反对称的周期延拓是圆周共轭反对称序列。

4.循环卷积与线性卷积
a. 计算序列x1[n]和x2[n]的线性卷积y[n]
b. 计算x1[n] 和x2[n]的傅里叶变换X1[k]和X2[k]，Y[k]=X1[k]*X2[k]，求Y[k]的反傅
里叶变换y2[n]（y2[n]即为x1[n]和x2[n]的循环卷积），比较y[n]与y2[n]的异同.
程序清单： clc;
clear;
x1=[1 3 5 3 6 8 3 9];
x2=[2 4 3 6 7 9 0 2 ];
y1=conv(x1,x2);
X1=fft(x1) X2=fft(x2) Y=X1.*X2 y2=ifft(Y) subplot(211);stem(y1);title('y1'); subplot(212);stem(y2);title('y2');
图像：
用Matlab 计算得:
X1 =-5.7071 - 9.7782i -1.0000 - 1.0000i -4.2929 - 5.7782i
X2 = 1.3640 - 3.7071i 6.0000 + 5.0000i -11.3640 + 2.2929i
Y = -0.0440 + 0.0078i -0.0010 - 0.0110i 0.0620 + 0.0558i
y2 = 170 175 178 143 161 115 154 158 4. 补零
用MATLAB 计算如下N 点序列的M 点DFT ： ()1010
n N x n ≤≤-⎧=⎨⎩其它
（1）取N=8，M=8
（2）取N=8，M=16
（3）取N=8，M=32
根据实验结果，分析同一序列不同点数的离散傅里叶变换的特点。

程序清单：
clc;
clear all;
close all;
N=8;
M=[8 16 32];
x=ones(1,N);
k_axis=[0:max(M)-1]; fft_x8=fft(x,M(1));
fft_x8=fftshift(fft_x8); fft_x16=fft(x,M(2));
fft_x16=fftshift(fft_x16);fft_x32=fft(x,M(3));
fft_x32=fftshift(fft_x32);
figure(1);
subplot(3,1,1);stem(k_axis(1:M(1)),abs(fft_x 8));title('8点序列的8点DFT');
subplot(3,1,2);stem(k_axis(1:M(2)),abs(fft_x 16));title('8点序列的16点DFT');
subplot(3,1,3);stem(k_axis(1:M(3)),abs(fft_x 32));title('8点序列的32点DFT');
图像：
结论：点数越多，傅里叶变换之后谱线间的间隔越小，也就是说分辨率越好。

实验二 时域采样与频域采样
1. 时域采样定理的验证，分析采样序列的特性
对x a (t )=A e -at sin(Ω0t)u (t )进行采样，可得到采样序列
x a (n)=x a (nT)=A e -anT sin(Ω0nT )u (n ),
其中A 为幅度因子，a 为衰减因子，Ω0是模拟角频率，T 为采样间隔。

a. 取采样频率f s =1 k Hz, 即T =1 ms ，观察|X(e j ω)|。

b. 改变采样频率，f s =300 Hz ，观察|X(e j ω)|的变化；
c.进一步降低采样频率，f s =200 Hz ，观察频谱混叠是否明显存在， 说明原因，
并记录|X(e j ω)|曲线。

程序清单： n=0:63; title('f=1000Hz ，|X(ej ω)|');
044.128,502,502==Ω=A a ππ
f1=1000;
f2=300;
f3=200; T1=1/f1;
T2=1/f2;
T3=1/f3;
a=50*sqrt(2)*pi;
xa=444.128*exp(-a*n*T1).*sin(a *n*T1);
Xa=T1*fft(xa,64);
subplot(321);stem(abs(xa)); title('f=1000Hz ，xa(n)'); subplot(322);plot(abs(Xa)); xb=444.128*exp(-a*n*T2).*sin(a*n*T2); Xb=T2*fft(xb,64); subplot(323);stem(abs(xb)); title('f=300Hz ，xa(n)'); subplot(324);plot(abs(Xb)); title('f=300Hz ，|X(ej ω)|'); xc=444.128*exp(-a*n*T3).*sin(a*n*T3); Xc=T3*fft(xb,64); subplot(325);stem(abs(xc)); title('f=200Hz ，xa(n)'); subplot(326);plot(abs(Xc)); title('f=200Hz ，|X(ej ω)|'); 图像：
2. 频域采样定理的验证。

给定信号
1013()2714260others n n x n n n +≤≤⎧⎪=-≤≤⎨⎪⎩
编写程序分别对频谱函数()[]()j X e DTFT x n ω=在区间[]0,2π上等间隔采样32点和16点，得到X 32(k)和X 16(k)，再分别对其进行32点和16点的IFFT ，得到x 32(n)和x 16(n)，分别画出()j X e ω、X 32(k)和X 16(k)的幅度谱。

并绘图显示x(n)、x 32(n)和x 16(n)的图形，进行对比分析，验证频域采样定理。

程序清单：
M=27;
N=32; n=0:M; k=0:N/2-1;
xa=(0:floor(M/2))+1; xb=ceil(M/2)-1:-1:0; xn=[xa,xb];
Xn=fft(xn,2048); Xk32=fft(xn,32); xn32=ifft(Xk32); X16k=Xk32(1:2:N);
x16n=ifft(X16k,N/2); Ab1=abs(Xn); Ab32=abs(Xk32); Ab16=abs(Xk16);
subplot(3,2,1),stem(Ab1);title('a'); subplot(3,2,2),stem(xn);title('a'); subplot(3,2,3),stem(Ab32);title('a'); subplot(3,2,4),stem(xn32);title('a'); subplot(3,2,5),stem(Ab16);title('a'); subplot(3,2,6),stem(x16n);title('a');
图像：
实验三 IIR 滤波器设计
一、
实验内容
采样频率为1Hz ，设计一个Butterworth 低通数字滤波器，其中通带临界频率0.2p f Hz =，通带内衰减小于1dB(1p dB α=)，阻带临界频率0.3s f Hz =，阻带内衰减大于25dB(25s dB α=)。

求这个数字滤波器的传递函数H(z)，输出它的幅频特性曲线。

利用冲激响应不变法和双线性变换法实现该滤波器，并将结果进行比较。

程序清单：
冲激响应不变法 Wp=2*pi*0.2; Ws=2*pi*0.3; Rp=1; Rs=25;
f1=linspace(0,Wp,5); f2=linspace(Wp,Ws,15); f3=linspace(Ws,2*pi*10,30);
T=1; Fs=1/T;
omegap=(2/T)*tan(Wp/2);
omegas=(2/T)*tan(Ws/2);
[N,Wn]=buttord(omegap,omegas,Rp,Rs,'s');
[B,A]=butter(N,Wn,'s');
h1=20*log10(abs(freqs(B,A,f1))); h2=20*log10(abs(freqs(B,A,f2))); h3=20*log10(abs(freqs(B,A,f3))); plot([f1f2 f3]/(2*pi),[h1,h2,h3]); grid;
xlabel('Frequency in Hz'); ylabel('Gain in dB');
图像：
程序清单：
双线性变换法设计 clc; clear; Wp=2*pi*0.2; Ws=2*pi*0.3;
[N,Wn]=buttord(omegap,omegas,Rp,Rs,'s'); [B,A]=butter(N,Wn,'s'); [b,a]=bilinear(B,A,Fs); [h,w]=freqz(b,a,256);
Rp=1; Rs=25; T=1;
Fs=1/T;
omegap=(2/T)*tan(Wp/2); omegas=(2/T)*tan(Ws/2);
h1=20*log10(abs(h)); plot(w/pi,h1); grid;
xlabel('Frequency in Hz'); ylabel('Gain in dB');
图像：
实验四 FIR 滤波器设计
实验内容：
用窗函数法设计一个长度N=8的线性相位FIR 滤波器。

其理想的幅频特性为
()100.400.4ωωπ
πωπ≤<⎧=⎨≤≤⎩
j d H e
分别用矩形窗、Hanning 窗、Hamming 窗、Blackman 窗设计该滤波器，并比较设计结果；
如果N=15，重复这一设计，观察幅频特性和相位特性的变化，注意长度N 变化对结果的影响。

程序清单：
用窗函数法设计一个长度N=8的线性相位FIR 滤波器。

clc;
Window=boxcar(8);
b=fir1(7,0.4,Window)
freqz(b,1);
图像
用窗函数法设计一个长度N=15的线性相位FIR 滤波器。

clc;
Window=boxcar(15); b=fir1(14,0.4,Window); freqz(b,1);
用Hanning 窗法设计一个长度N=8的线性相位FIR 滤波器。

clc;
Window=hanning(8);
b=fir1(7,0.4,Window); freqz(b,1);
用Hanning 窗法设计一个长度N=15的线性相位FIR 滤波器。

clc;
b=fir1(14,0.4,Window);
Window=hanning(15); freqz(b,1);
用Hamming窗法设计一个长度N=8的线性相位FIR滤波器。

clc;
Window=hamming(8); b=fir1(7,0.4,Window); freqz(b,1);
用Hamming窗法设计一个长度N=15的线性相位FIR滤波器。

clc;
Window=hamming(15); b=fir1(14,0.4,Window); freqz(b,1);
用Blackman 窗法设计一个长度N=8的滤波器。

clc;
Window=blackman(8);
b=fir1(7,0.4,Window); freqz(b,1);
用Blackman 窗法设计一个长度N=15的滤波器。

clc;
Window=blackman(15);
b=fir1(14,0.4,Window); freqz(b,1);
2.人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经多低通滤波处理后磁能作为判断心脏功能的有用信息。

下面给出一实际心电图信号采样序列样本x (n )，其中存在高频干扰，采用上述设计的某个低通滤波器滤除其中的高频干扰成分，在同一幅图中绘出滤波前后的序列x (n ) 和y (n )。

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,
0,-16,-38,-60,-84,-90,-66,-32,-4 ,-2,-4,8,12, 12,10, 6, 6, 6, 4, 0, 0,0,0, 0, -2 ,-4,0,0, 0, -2, -2,0,0,-2,-2,-2,-2,0];
Window=boxcar(8); b=fir1(7,0.4,Window);
y=filter(b,1,x);
plot(x);hold on;plot(y,'r');
图像：。