数字信号处理实验集合

实验一 信号、系统及系统响应
一、实验目的
1、 熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。
2、 熟悉时域离散系统的时域特性。
3、 利用卷积方法观察分析系统的时域特性。
4、 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验代码

%dsp11

clear all;
close all;
figure(1)
a1 = FF(444.128,222,222,1000);
n = 0:49;
[b1,K,D] = dtft(a1,n,50);
subplot(3,2,1);%3行2列 1 位置
stem([0:49],a1,'.')%[]:n的范围  a1：纵坐标  画  离散图
xlabel('n')
ylabel('xa(n)');
title('fs=1000');
subplot(3,2,2)
plot(K*D,abs(b1))%： 连续图
xlabel('w')
ylabel('|X(jw)|');
title('fs=1000');

a2 = FF(444.128,222,222,300);
n = 0:49;
[b2,K,D] = dtft(a2,n,50);
subplot(3,2,3);%3行2列 1 位置
stem([0:49],a1,'.')%[]:n的范围  a1：纵坐标  画  离散图
xlabel('n')
ylabel('xa(n)');
title('fs=300');
subplot(3,2,4)
plot(K*D,abs(b2))%： 连续图
xlabel('w')
ylabel('|X(jw)|');
title('fs=300');

a3 = FF(444.128,222,222,200);
n = 0:49;
[b3,K,D] = dtft(a3,n,50);
subplot(3,2,5);%3行2列 1 位置
stem([0:49],a1,'.')%[]:n的范围  a1：纵坐标  画  离散图
xlabel('n')
ylabel('xa(n)');
title('fs=200');
subplot(3,2,6)
plot(K*D,abs(b3))%： 连续图
xlabel('w')
ylabel('|X(jw)|');
title('fs=200');

%dsp12

clear all;
close all;
figure(2)
[Xb,n1]=impseq(0,0,0);
[Ha,n2]=stepseq(0,0,9);
Hb=impseq(0,0,3)+2.5*impseq(1,0,3)+2.5*impseq(2,0,3)+impseq(3,0,3);

subplot(3,2,1);
stem(n1,Xb,'.');
xlabel('n');
ylabel('xb(n)');
N=64;
[A1,k1,D1]=dtft(Xb,n1,N);
subplot(3,2,2),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Xb(jw)|');

subplot(3,2,3);
stem([0:3],Hb,'.');
xlabel('n');
ylabel('hb(n)');
N=64;
[A1,k1,D1]=dtft(Hb,0:3,N);
subplot(3,2,4),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Hb(jw)|');

Xb=double(Xb);
Hb=double(Hb);
y=conv(Hb,Xb);
subplot(3,2,5);
stem([0:3],y,'.');
xlabel('w');
ylabel('yb(n)');

N=64;
[A1,k1,D1]=dtft(y,0:3,N);
subplot(3,2,6),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Y1(jw)|');

figure(3)
Ha=double(Ha);
y=conv(Ha,Ha);
subplot(2,2,1);
stem(n2,Ha,'.');
xlabel('n');
ylabel('Ha(n)');

N=64;
[A1,k1,D1]=dtft(Ha,n2,N);
subplot(2,2,2),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Ha(jw)|');

subplot(2,2,3);
stem([0:18],y,'.');
xlabel('n');
ylabel('y(n)');

N=64;
[A1,k1,D1]=dtft(y,0:18,N);
subplot(2,2,4),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Y2(jw)|');

figure(4)
[c,n4]=stepseq(0,0,4);
Ha=double(Ha);
y=conv(Ha,c);
subplot(2,2,1);
stem(n4,c,'.');
xlabel('n');
ylabel('Ha(n)');

N=64;
[A1,k1,D1]=dtft(c,n4,N);
subplot(2,2,2),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Xc(jw)|');

subplot(2,2,3);
stem([0:13],y,'.');
xlabel('n');
ylabel('y(n)');

N=64;
[A1,k1,D1]=dtft(y,0:13,N);
subplot(2,2,4),plot(k1*D1),abs(A1);
xlabel('w');
ylabel('|Y(jw)|');

%dsp13

clear all;
close all;
figure(5)
Xa = FF(4, 0.4, 2.073, 1);
Xb = impseq(0, 0, 1);
Hb = impseq(0, 0, 3) + 2.5 * impseq(1, 0, 3) + 2.5 * impseq(2, 0, 3) + impseq(3, 0, 3);
b1 = abs(fft(Xa));
subplot(3, 2, 1);
stem([0:49], Xa, '.');
xlabel('n');
ylabel('xa(n)');
title('fs = 1');

N = 512;
[A1, k1, D1] = dft2(Xa, N);
subplot(3, 2, 2), plot(k1 * D1, abs(A1));
xlabel('w');
ylabel('|Xa(jw)|');

subplot(3, 2, 3);
stem([0:3], Hb, '.');
xlabel('n');
ylabel('Hb(n)');

N = 512;
[B1, k1, D1] = dft2(Hb, N);
subplot(3, 2, 4), plot(k1 * D1, abs(B1));
xlabel('w');
ylabel('|Hb(jw)|');

Hb = double(Hb);
Xa = double(Xa);
y = conv(Hb, Xa);
%Y = abs(fft(y, 64));
subplot(3, 2, 5);
stem([0:52], y, '.');
xlabel('n');
ylabel('y(n)');

N = 512;
[Y1, k1, D1] = dft2(y, N);
subplot(3, 2, 6), plot(k1 * D1, abs(Y1));
xlabel('w');
ylabel('|Y(jw)|');

figure(6)
XH = A1 .* B1;
yi = ifft(XH, 512);

subplot(2, 1, 1);
stem([0:52], y, '.');
xlabel('n');
ylabel('y(n)');

subplot(2, 1, 2);
stem([0:52], real(yi(1:53)), '.');
xlabel('n');
ylabel('yi(n)');

%impseq
function [x, n] = impseq(n0, n1, n2)
n = [n1:n2];
x = [(n - n0) == 0];

%stepseq
function [x,n] = stepseq(n0,n1,n2)
n = [n1:n2];
x = [(n-n0) >= 0];

%FF
function c = FF(A,a,w,fs)
n = 0: 49;
c = A*exp((-a)*n/fs).*sin(w*n/fs).*stepseq(0,0,49);

%dtft
function [X,K,D] = dtft(x,n,N)%K:离散频谱的点的序号 D ：离散频谱点的间隔 X：离散频谱
K = floor((-N/2+0.5):(N/2-0.5));
D = 2*pi/N;
X= x*exp(-j*D*n'*K);%n’转之

实验运行结果截图

实验二 用FFT作谱分析
一、实验目的
1、 进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。
2、 熟悉FFT算法原理和FFT子程序的应用。
3、 学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
二、实验代码

%dsp2
clear all;
close all;
x1=[1,1,1,1];
x2=[1,2,3,4,4,3,2,1];
x3=[4,3,2,1,1,2,3,4];
x6=xinhao6(16,32);
N1=8;
n=0:N1-1;
x4=cos(0.25*pi*n);
N2=16;
n=0:N2-1;
x5=sin(0.125*pi*n);

figure(1);
[X1a,K1,D1]= dft1(x1,8);
[X1b,K2,D2]=dft1(x1,64);

subplot(3,1,1);
stem([0:3],x1,'.')
ylabel('x1(n)');

subplot(3,1,2);
stem(K1*D1,abs(X1a),'.');
xlabel('w');
ylabel('|X1a|(N=8)');

subplot(3,1,3);
stem(K2*D2,abs(X1b),'.');
xlabel('w');
ylabel('|X1a|(N=8)');

figure(2);
[X2a,K1,D1]= dft1(x2,8);
[X2b,K2,D2]=dft1(x2,16);

subplot(3,1,1);
stem([0:7],x2,'.')
ylabel('x2(n)');

subplot(3,1,2);
stem(K1*D1,abs(X2a),'.');
xlabel('w');
ylabel('|X2a|(N=8)');

subplot(3,1,3);
stem(K2*D2,abs(X2b),'.');
xlabel('w');
ylabel('|X2b|(N=16)');

figure(3);
[X3a,K1,D1]= dft1(x3,8);
[X3b,K2,D2]=dft1(x3,16);

subplot(3,1,1);
stem([0:7],x3,'.')
ylabel('x3(n)');

subplot(3,1,2);
stem(K1*D1,abs(X3a),'.');
xlabel('w');
ylabel('|X3a|(N=8)');

subplot(3,1,3);
stem(K2*D2,abs(X3b),'.');
xlabel('w');
ylabel('|X3b|(N=16)');

figure(4);
[X4,K1,D1]= dft1(x4,8);

subplot(3,1,1);
stem([0:7],x4,'.')
ylabel('x4(n)');

subplot(3,1,2);
stem(K1*D1,abs(X4),'.');
xlabel('w');
ylabel('|X4|(N=8)');

figure(5);
[X5,K1,D1]= dft1(x5,16);

subplot(3,1,1);
stem([0:15],x5,'.')
ylabel('x5(n)');

subplot(3,1,2);
stem(K1*D1,abs(X5),'.');
xlabel('w');
ylabel('|X5|(N=16)');

figure(6);
[X6,K1,D1]= dft1(x6,16);

subplot(3,1,1);
stem([0:15],x6,'.')
ylabel('x6(n)');

subplot(3,1,2);
stem(K1*D1,abs(X6),'.');
xlabel('w');
ylabel('|X6|(N=16)');

N3=8;
n=0:N3-1;
x4=cos(0.25*pi*n);
x5=sin(0.125*pi*n);
x7=x4+x5;
x8=x4+x5*i;

figure(7);
[X7a,K1,D1]= dft1(x7,8);
[X7b,K2,D2]=dft1(x7,16);

subplot(3,1,1);
stem([0:7],x7,'.')
ylabel('x7(n)');

subplot(3,1,2);
stem(K1*D1,abs(X7a),'.');
xlabel('w');
ylabel('|X7a|(N=8)');

subplot(3,1,3);
stem(K2*D2,abs(X7b),'.');
xlabel('w');
ylabel('|X7b|(N=16)');

figure(8);
[X8a,K1,D1]= dft1(x8,8);
[X8b,K2,D2]=dft1(x8,16);

subplot(4,1,1);
stem([0:7],real(x8),'.')
ylabel('x8(n)实部');

subplot(4,1,2);
stem(K1*D1,imag(x8),'.');
ylabel('x8(n)虚部');

subplot(4,1,3);
stem(K1*D1,abs(X8a),'.');
xlabel('w');
ylabel('|X8a|(N=8)');

subplot(4,1,4);
stem(K2*D2,abs(X8b),'.');
xlabel('w');
ylabel('|X8b|(N=16)');

%figure(9）；
%subplot（2，1，1）；
%stem(K1*D1,real(X8),'.');
%xlabel('w');
%ylabel('X8(real)(N=8)');
%
%subplot(2,1,2);
%stem(K1*D1,imag(X8),'.');
%xlabel('w');
%ylabel('X8(image)(N=8)');
%rr=real(X8);
%ii=imag(X8);

%dft1
function [X,K,D]=dft1(x,N)
X=fftshift(fft(x,N));
D=2*pi/N;
K=floor((-N/2+0.5):(N/2-0.5));

%xinhao6
function x=xinhao6(N,fs)
n=0:N-1;
x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);

实验运行结果截图


实验三 用双线性变换法设计IIR数字滤波器

一、实验目的
1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。
2、掌握数字滤波器的计算机仿真方法。
3、通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。
二、实验代码

%dsp3_a
clear all;
close all;
%数字指标
Wp = 0.2 * pi;%通带截止频率
Ws = 0.3 * pi;%阻带截止频率
Rp = 1;%通带最大衰减单位dB
Rs = 15;%阻带最小衰减单位dB
%转换成模拟域指标

T = 1;Fs = 1 / T;
Wap = (2 / T )* tan(Wp / 2);
Was = (2 / T ) * tan(Ws / 2);

%模拟滤波器设计
[n, Wn] = buttord(Wap, Was, Rp, Rs, 's');%获得模拟滤波器阶数和通带截止频率
[z, p, k] = buttap(n);%设计模拟归一化原型低通滤波器，求得其零，极点和增益
[b0, a0] = zp2tf(z, p, k);%将零，极点转化为模拟滤波器的传递函数表示
[b, a] = lp2lp(b0, a0, Wn);%将模拟归一化原型低通滤波器转化为要求设计的去归化一的模拟低通滤波器

%使用双线性变换法得到设计的数字滤波器
[bz, az] = bilinear(b, a, Fs);


[h, w] = freqz(bz, az);%计算频率响应

%绘幅频率响应曲线
figure(1)
plot(w, abs(h));
figure(2)
plot(w, 20 * log10(abs(h)));

%信号滤波
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];
figure(3);
subplot(2, 1, 1);
stem(x);
y = filter(bz, az, x);
subplot(2, 1, 2);
stem(y);

%dsp3_b
clear all;
close all;

Wp = 0.2;
Ws = 0.3;
Rp = 1;
Rs = 15;

[n,Wn] = buttord(Wp,Ws,Rp,Rs);  %计算模拟滤波器阶数N和频率
[b,a]= butter(n,Wn);
[h,w] = freqz(b,a);   %计算频率响应

figure(1); %幅频特性曲线
plot(w,abs(h));
axis([0 pi/2 0 1]);%设置x，y轴显示界限
figure(2);  %幅频特性曲线（dB值）
plot(w,20*log10(abs(h)));

%信号滤波
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];
figure(3);
subplot(2, 1, 1);
stem(x);
y = filter(b, a, x);
subplot(2, 1, 2);
stem(y);

实验运行结果截图

实验四 用窗函数法设计FIR数字滤波器

一、实验目的：
1、掌握用窗函数法设计FIR数字滤波器的原理和方法。
2、熟悉线性相位FIR数字滤波器特性。
3、了解各种窗函数对滤波特性的影响。
二、实验代码：

%dsp4_a
clear all;
close all;
%升余弦窗设计线性低通FIR数字滤波器，截止频率wc=pi/4，窗口长度N=15,33
wc=0.25*pi;N1=15;N2=33;
hd1=ideallp(wc,N1);
hd2=ideallp(wc,N2);
wd1=hanning(N1)'; %'表示转置
b1=hd1.*wd1;
wd2=hanning(N2)'; 
b2=hd2.*wd2;
[H1,w]=freqz(b1,1);
[H2,w]=freqz(b2,1);

figure(1)
subplot(2,1,1)
stem([0:N1-1],b1,'.');
title('h1(n)');
subplot(2,1,2)
stem([0:N2-1],b2,'.');
title('h2(n)');

%幅频响应
figure(2)
%一个窗口两条曲线，b代表blue，r代表红色(也可用不同线型画曲线)
plot(w,abs(H1),'b',w,abs(H2),'r');
xlabel('w');
ylabel('|H(w)|');
legend('升余弦窗N=15','升余弦窗N=33');%加图例标识，右上方呈现

%对数幅频特性曲线图
figure(3)
plot(w,20*log10(abs(H1)),'b',w,20*log10(abs(H2)),'r');
xlabel('w');
ylabel('dB');
legend('升余弦窗N=15','升余弦窗N=33');

%相频响应
figure(4)
plot(w,phase(H1),'b',w,phase(H2),'r');%phase：对频率响应取相位响应
xlabel('w');
ylabel('相位响应');
legend('升余弦窗N=15','升余弦窗N=33');

%dsp4_b
clear all;
close all;
%用矩形窗、升余弦窗、改进升余弦窗、二阶升余弦窗设计线性相位低通滤波器
wc=0.25*pi;
N=33;
hd=ideallp(wc,N);
wd1=boxcar(N)';   b1=hd.*wd1;   %矩形窗
wd2=hanning(N)';  b2=hd.*wd2;   %升余弦窗
wd3=hamming(N)';  b3=hd.*wd3;   %改进升余弦窗
wd4=blackman(N)'; b4=hd.*wd4;   %二阶升余弦窗
[H1,w]=freqz(b1,1);
[H2,w]=freqz(b2,1);
[H3,w]=freqz(b3,1);
[H4,w]=freqz(b4,1);

figure(1)
plot([0:N-1],wd1,'r',[0:N-1],wd2,'y',[0:N-1],wd3,'b',[0:N-1],wd4,'g');
legend('矩形窗','升余弦窗','改进升余弦窗','二阶升余弦窗');
title('w(n)');

figure(2)
subplot(4,1,1)
stem([0:N-1],b1,'.');
title('h1(n)');
subplot(4,1,2)
stem([0:N-1],b2,'.');
title('h2(n)');
subplot(4,1,3)
stem([0:N-1],b3,'.');
title('h3(n)');
subplot(4,1,4)
stem([0:N-1],b2,'.');
title('h4(n)');

%幅频响应
figure(3)
plot(w,abs(H1),'r',w,abs(H2),'y',w,abs(H3),'b',w,abs(H4),'g');
xlabel('w');
ylabel('|H(w)|');
legend('矩形窗','升余弦窗','改进升余弦窗','二阶升余弦窗');

%对数幅频响应
figure(4)
plot(w,20*log10(abs(H1)),'r',w,20*log10(abs(H2)),'y',w,20*log10(abs(H3)),'b',w,20*log10(abs(H4)),'g');
xlabel('w');
ylabel('dB');
legend('矩形窗','升余弦窗','改进升余弦窗','二阶升余弦窗');

%相频响应
figure(5)
plot(w,phase(H1),'r',w,phase(H2),'y',w,phase(H3),'b',w,phase(H4),'g');
xlabel('w');
ylabel('相位响应');
legend('矩形窗','升余弦窗','改进升余弦窗','二阶升余弦窗');

%ideallp
%求理想低通滤波器的单位脉冲响应hd
function hd=ideallp(wc,N) %wc=截止频率，N=滤波器长度
tao=(N-1)/2;
n=[0:(N-1)];
m=n-tao+eps;   %加上一个极小数（eps）以避免0做除数
hd=sin(wc*m)./(pi*m);

实验运行结果截图