《数字信号处理》Matlab实验
一．离散信号的 FFT 分析
知识点：利用FFT对信号频谱进行分析，用DFT进行信号分析时基本参数的选择，以及信号经过离散时间傅立叶变换（DTFT）和有限长度离散傅立叶变换（DFT）后信号频谱上的区别。
实验教学内容：
1.用Matlab编程上机练习。已知：
N=25。这里Q=0.9+j0.3。可以推导出 ，
首先根据这个式子计算X(k)的理论值，然后计算输入序列x(n)的32个值，再利用基2时间抽选的FFT算法，计算x(n)的DFT X(k)，与X(k)的理论值比较（要求计算结果最少6位有效数字）。
解:
format long
Q=0.9+0.3i;
WN=exp(-2*pi*1i/32);
Xk=(1-Q^32)./(1-Q*WN.^[0:24]);
xn=Q.^[0:24];
Xkfft=fft(xn,32);
for (k0=1:1:25)
difference=Xk(k0)-Xkfft(k0);
end;
subplot(3,1,1);stem(abs(Xk(1:1:24)),'.');title('DFT x(n)');xlabel('k');axis([0,35,0,15]);
subplot(3,1,2);stem(abs(Xkfft(1:1:32)),'g.');title('FFT x(n)');xlabel('k');axis([0,35,0,15]);
subplot(3,1,3);stem(abs(difference(1:1:25)),'r.');title('Xk-Xkfft');xlabel('k');axis([0,35,0,15]);
2. 假设信号 x(n) 由下述信号组成：
请选择合适的长度 N 和窗函数，用 DFT 分析其频谱，得到清楚的三根谱线。
解：
N = 1000; % Length of DFT
n = [0:1:N-1];
xn = 0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-pi/4);
Xk = fft(xn,N);
xn1=0.001*cos(0.45*n*pi);
Xk1 = fft(xn1,N);
xn2=sin(0.3*n*pi);
Xk2 = fft(xn2,N);
xn3=-cos(0.302*n*pi-pi/4);
Xk3 = fft(xn3,N);
k=[0:1:N/2];
subplot(5,1,1);stem(k,abs(xn(1:1:(N/2+1))),'.');title('x(n)');xlabel('k');axis([140,240,0,4]);
subplot(5,1,2);stem(k,abs(Xk(1:1:(N/2+1))),'r.');title('DFT 0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-pi/4)');xlabel('k');axis([140,240,0,4]);
subplot(5,1,3);stem(k,abs(Xk1(1:1:(N/2+1))),'g.');title('DFT 0.001*cos(0.45*n*pi)');xlabel('k');axis([140,240,0,4]);
subplot(5,1,4);stem(k,abs(Xk2(1:1:(N/2+1))),'b');title('DFT sin(0.3*n*pi)');xlabel('k');axis([140,240,0,4]);
subplot(5,1,5);stem(k,abs(Xk3(1:1:(N/2+1))));title('DFT -cos(0.302*n*pi-pi/4)');xlabel('k');axis([140,240,0,4]);
二．IIR 数字滤波器的设计和实现
DTMF 信号的编码: 把你的联系电话号码 DTMF 编码生成为一个 .wav 文件。
技术指标：
根据 ITU Q.23 建议，DTMF 信号的技术指标是：传送/接收率为每秒 10 个号码，或每个号码 100ms。每个号码传送过程中，信号存在时间至少 45ms，且不多于 55ms，100ms 的其余时间是静音。在每个频率点上允许有不超过 ±1.5% 的频率误差。任何超过给定频率 ±3.5% 的信号，均被认为是无效的，拒绝接收。
代码如下：
clc
for(K=1:10)
clear
phone=input('请输入电话号码','s');
sum=length(phone);
fs=8000;
total_x=[];
fH=[1336,1209,1336,1477,1209,1336,1477,1209,1336,1477];
fL=[941,697,697,697,770,770,770,852,852,852];
for i=1:sum
num=phone(i)-48+1;
n=1:400;
x=sin(2*pi*n*fL(num)/fs)+ sin(2*pi*n*fH(num)/fs);
x=[x,zeros(1,400)];
total_x=[total_x,x];
end
wavwrite(total_x,'num_sound')
sound(total_x);
plot(total_x);axis([0,10000,-2,2]);
title('电话号码时域波形')
k=[18,20,22,25,32,35,38];
N=210;
tm=[49,50,51;52,53,54;55,56,57;0,48,0];
for i=1:sum
j=800*(i-1);
X=goertzel(total_x(j+1:j+N),k+1);
value=abs(X);
figure(2)
subplot(2,6,i);
stem(k,value,'.','r');
title('FFT x(n)');xlabel('k');ylabel('|X(k)|');
limit=20;
for i1=5:7
if value(i1)>limit
break;
end
end
for j1=1:4
if value(j1)>limit
break;
end
end
buffer(i)=tm(j1,i1-4);
end
disp(['接收端检测到的号码'])
disp(setstr(buffer))
end
三．FIR 数字滤波器的设计和实现
录制自己的一段声音，长度为 45秒，取样频率 32kHz，然后叠加一个高斯白噪声，使得信噪比为 20dB。请采用窗口法设计一个 FIR 带通滤波器，滤波噪声提高质量。请随实验报告同时提交原始加噪和滤波后的声音文件（.mp3 格式）。
提示：
1. 滤波器指标参考：通带边缘频率为 4kHz，阻带边缘频率为4.5kHz，阻带衰减大于 50dB；
2. Matlab 函数 y = awgn(x,snr,'measured') ，首先测量输入信号 x 的功率，然后对其叠加高斯白噪声；
第二篇：北邮《数字信号处理》Matlab实验——周期序列的谱分析
实验题目2：周期序列的谱分析?实验目的：?
利用DFT分析模拟信号xa?t??cos?16?t?之频谱。?
?
试验内容：?
1、设定采样周期T并说明原因?
根据奈奎斯特取样条件，fs>=2f=2*8=16Hz，所以取fs=20Hz，T=1/f=0.05s。??
2、若令x(n)=cos(16πnT)，确定该序列之周期N并说明原因?因为x(n)?cos(42?5?，所以周期N=5. ?n)，5?2
5
?
3、绘制10个周期内x(n)的取值情况?
T = 1/20; t = 0:T:5-T; N = 5;
figure(1);
xn = cos(16*pi*t); stem(0:10*N-1, xn(1:10*N));
title('10个周期内x(n)的取值情况'); ylabel('x(n)'); xlabel('n');
?
4、令x1(n)表示x(n)的主值序列，绘制|DFT(x1(n)|，解释取值情况?figure(2);
x1n = xn(1:N); stem(0:N-1, x1n);
title('x(n)的主值序列'); ylabel('x1(n)'); xlabel('n'); figure(3);
Xk = fft(x1n, N); stem(0:N-1, abs(Xk)); title('|DFT(x1(n))|'); ylabel('X1(k)'); xlabel('k');
?
?
所得X1(k)取值即为x(n)的5点（一个周期内）DFT结果，为清晰谱线。??
5、令x2(n)表示x(n)的任一周期，绘制|DFT(x2(n)|，解释取值情况?figure(4);
r = 4; x2n = xn(r:r+N-1); stem(0:N-1, x2n);
title('x(n)的任一周期'); ylabel('x2(n)'); xlabel('n'); figure(5);
Xk = fft(x2n, N); stem(0:N-1, abs(Xk));
title('|DFT(x2(n))|'); ylabel('X2(k)'); xlabel('k');
?
?
所得X2(k)取值与x(n)主值序列DFT结果相同。因为DFT是反映信号的频域特性，所以同为一个周期，频域特性一定相同，无论起始位置如何，其DFT情况都相同。??
6、令x3(n)表示x(n)的2个周期，绘制|DFT(x3(n)|，解释取值情况?
figure(6);
x3n = xn(1:2*N); stem(0:2*N-1, x3n);
title('x(n)的2个周期'); ylabel('x3(n)'); xlabel('n'); figure(7);
Xk = fft(x3n, 2*N); stem(0:2*N-1, abs(Xk)); title('|DFT(x3(n))|'); ylabel('X3(k)'); xlabel('k');
?
?
所得X3(k)是在x(n)主值序列DFT结果的基础上插入取值为零的点，将点数提高至10，仍为清晰谱线。因为选取了2个周期，所以DFT点数是X1(k)的2倍。因为取的周期的整数倍，所以可以得到清晰谱线。?
?
7、若x4(n)=x(n)×RM(n)，而M不是x(n)周期的整数倍，绘制|DFT(x4(n)|，解释取值情况?
figure(8);
M = 8; x4n = xn(1:M); stem(0:M-1, x4n);
title('x(n)*R8(n)'); ylabel('x4(n)'); xlabel('n');
figure(9);
Xk = fft(x4n, M); stem(0:M-1, abs(Xk));
title('|DFT(x4(n))|'); ylabel('X4(k)'); xlabel('k');
?
?
所得X4(k)是非周期整数倍点数的DFT，不是清晰谱线。因为选取了1.3个周期，不是周期的整数倍，选取的x4(n)没有体现出x(n)的周期特性，所以得到的谱线的形状特征必然与X1(k)、X2(k)、X3(k)不同。?
?