数字信号处理2： 离散信号与系统的频谱分析
文章目录
- 前言
-
第一部分 实验目标
-
第二部分 实验设备
-
第三部分 实验内容
-
第四部分 实验原理
-
第五部分 实验步骤
- 1.基于序列的离散傅里叶变换分析
- 2.利用共轭对称性设计高效算法来计算两个N点实序列的DFT
- 3.研究线性卷积与循环卷积的实现及其相互关系
- 4.对比DFT与FFT运算时间效率的研究
- 5.探讨FFT在信号频谱分析中的应用及其对采样频率和噪声的影响研究
- 6.拓展创新训练内容
- 6.1研究信号持续时间、频谱分析范围、采样点数与谱分辨率之间的关系研究
- 6.2探讨频域采样恢复技术在频谱内插中的应用探讨
- 6.3完成语音信号的基本特性分析任务。
-
六、实验结论与心得体会（手写）
-
七、实验参考资料
-
前言
本人初次学习离散信号与系统的频谱分析。
一、实验目的
- 深入理解离散傅里叶变换（DFT）及其快速算法（FFT）的具体实现方案。
- 分析序列DFT的基本特性。
- 熟练掌握通过DFT/FFT计算序列线性卷积的具体方法。
- 深入研究如何通过DFT/FFT计算序列线性卷积的具体方法；全面了解可能出现的分析误差及其对结果的影响，并在此基础上正确应用这些方法。
- 探讨采样频率对谱分析结果的影响机制。
- 深入理解并掌握利用快速傅里叶变换（FFT）进行语音信号分析的方法。
二、实验设备
1.计算机
2.Matlab软件2014a以上版本。
三、实验内容
- 分别对待不同输入序列并对其执行离散傅里叶变换后进行深入分析；DFT所具有的共轭对称性质得以在单次N点快速傅里叶变换（FFT）中有效体现出来以同时计算两个独立的N点实数序列的DFT值。
- 线性卷积与循环卷积之间存在本质区别即前者是无限长序列而后者是有限长周期序列；同时可借助快速傅里叶变换算法来实现线性卷积过程从而避免直接计算所带来的高复杂度。
- 通过对比计算不同长度信号的DFT与FFT运算量评估它们之间的算法效率差异进而优化数据处理流程以提高整体性能指标。
- 基于快速傅里叶变换技术可有效实现带噪声信号特征提取过程从而降低信噪比提高检测灵敏度。
- 采用快速傅里叶变换方法能够精准地解析信号频谱从而为后续的数据分析提供可靠的基础依据。
- 本节内容主要探讨了离散系统中采样频率时域持续时间和频谱分辨率三者之间的相互影响关系此外还研究了基于内插技术恢复高分辨率频谱的方法并对语音信号进行了初步处理以揭示其固有特性
四、实验原理
- 信号的离散傅里叶变换及其特性
- 通过调用FFT算法完成这一操作其计算过程涉及调用FFT算法完成这一操作
- 通过调用FFT算法完成这一操作完成对信号频域特性的分析
- 前者是后者的约四分之一前者是后者的约四分之一
- 信噪比（S/N）表示的是输入信号与噪声之间的比例关系
五、实验步骤
1.序列的离散傅里叶变换及分析
依次处理复序列、实序列、实偶序列、实奇序列以及虚奇序列表现其离散傅里叶变换的结果特性，并对其实验现象展开深入探讨。实验中所选用的信号类型可由研究者自由选定
x1=ones(1,5);
x2=1:5;
x=x1+j*x2;
y1=fft(x,10);
y2=fft(x1,10);
subplot(1,2,1)
stem(0:9,abs(y1),'filled');
title('复序列DFT的幅度谱');
subplot(1,2,2)
stem(0:9,abs(y2),'filled');
title('实序列DFT的幅度谱');
c
结果：
B、设定实偶序列：x(n)= {5,4,3,2,1,0,1,2,3,4}，实奇序列：x(n)= {0,4,3,2,1,0,-1,-2,-3,-4}，虚奇序列：x(n)= j*{0,4,3,2,1,0,-1,-2,-3,-4};
程序：
x=[5,4,3,2,1,0,1,2,3,4];
y=fft(x,10);
subplot(3,3,1)
stem(0:9,abs(y),'filled');
title('实偶序列DFT的幅度谱');
subplot(3,3,2)
stem(0:9,real(y),'filled');
title('实偶序列DFT的实部');
subplot(3,3,3)
stem(0:9,imag(y),'filled');
title('实偶序列DFT的虚部');
x=[0,4,3,2,1,0,-1,-2,-3,-4];
y=fft(x,10);
subplot(3,3,4)
stem(0:9,abs(y),'filled');
title('实奇序列DFT的幅度谱');
subplot(3,3,5)
stem(0:9,real(y),'filled');
title('实奇序列DFT的实部');
subplot(3,3,6)
stem(0:9,abs(imag(y)),'filled');
title('实奇序列DFT的虚部');
x=j*[0,4,3,2,1,0,-1,-2,-3,-4];
y=fft(x,10);
subplot(3,3,7)
stem(0:9,y,'filled');
title('虚奇序列DFT的幅度谱');
subplot(3,3,8)
stem(0:9,real(y),'filled');
title('虚奇序列DFT的实部');
subplot(3,3,9)
stem(0:9,imag(y),'filled');0
title('虚奇序列DFT的虚部');
c
结果：
分析：
2.利用共轭对称性，设计高效算法计算2个N点实序列的DFT。
采用一个N点快速傅里叶变换（FFT）来计算两个长度均为N的实序列各自的N点离散傅里叶变换（DFT），并对比这两种方法的结果。
x1=[1,2,3,4,5,6];
x2=[1,2,1,3,1,4];
y1=fft(x1);
y2=fft(x2);
x=x1+j*x2;
y=fft(x,6);
[xec,xoc]=circevod(y);
subplot(2,2,1)
stem(0:5,abs(y1),'filled');
title('x_1(n)的直接六点DFT的幅度谱');
grid on;
subplot(2,2,2)
stem(0:5,abs(y2),'filled');
title('x_2(n)的直接六点DFT的幅度谱');
grid on;
subplot(2,2,3)
stem(0:5,abs(xec),'filled');
title('x_1(n)高效算法DFT的幅度谱');
grid on;
subplot(2,2,4)
stem(0:5,abs(xoc),'filled');
title('x_2(n)高效算法DFT的幅度谱');
grid on;
c
结果：
分析：
3.线性卷积及循环卷积的实现及二者关系分析
计算两个序列的线性卷积和循环卷积特性。其中循环卷积采用时域法和频域法两种方式来计算。设序列x₁长度为M点，序列x₂长度为N点，则循环卷积长度L分别取大于、等于以及小于(M + N - 1)值的情况进行研究分析。实验中选取自定义参数：设定数字序列x₁(n) = {1, 2, 3, 4} (0 ≤ n ≤ 5)，x₂(n) = {1, 2, 1, 3} (0 ≤ n ≤ 5)，并基于上述理论实现不同L值下的循环卷积运算，并与理论推导结果进行对比验证。
程序实现部分则采用Matlab编程环境完成相关算法开发。
figure(1)
x1=[1,2,3,4];
x2=[1,2,1,3];
subplot(3,2,1)
y1=circonv(x1,x2,4);
stem(0:3,abs(y1),'filled');
title('从时域计算4点循环卷积');
grid on;
subplot(3,2,3)
y1=circonv(x1,x2,7);
stem(0:6,abs(y1),'filled');
title('从时域计算7点循环卷积');
grid on;
subplot(3,2,5)
y1=circonv(x1,x2,12);
stem(0:11,abs(y1),'filled');
title('从时域计算12点循环卷积');
grid on;
subplot(3,2,2)
y1=fft(x1,4);
y2=fft(x2,4);
x=ifft(y1.*y2)
stem(0:3,abs(x),'filled');
title('从频域计算4点循环卷积');
grid on;
subplot(3,2,4)
y1=fft(x1,7);
y2=fft(x2,7);
x=ifft(y1.*y2)
stem(0:6,abs(x),'filled');
title('从频域计算7点循环卷积');
grid on;
subplot(3,2,6)
y1=fft(x1,12);
y2=fft(x2,12);
x=ifft(y1.*y2)
stem(0:11,abs(x),'filled');
title('从频域计算12点循环卷积');
grid on;
figure(2)
x1=[1,2,3,4];
x2=[1,2,1,3];
y=conv(x1,x2);
stem(0:6,y,'filled')
grid on;
title('线性卷积的结果');
c
结果：
分析：
4.比较DFT和FFT的运算时间
（1）自行决定生成特定类型的序列用于后续运算。
（2）选取一系列递增的DFT点数N值：64；256；1024；2048；依此类推。
（3）本实验将采用两种不同的方法进行DFT运算：一种是基于定义法的手动实现；另一种是调用FFT算法函数。借助tic-toc计时指令记录每次运算所需的时间，并详细分析两组数据之间的差异性；同时对比不同点数N下的运算时间变化及其效率提升情况。
解：构造测试数据集x为[1, 2, 3, 4]以及长度为60的一系列随机数。
程序：
x=[1,2,3,4,ones(1,60)];
tic
X1=dft(x,64);
toc
tic
X2=fft(x,64);
toc
时间已过 0.003188 秒。
时间已过 0.001967 秒。
x=[1,2,3,4,ones(1,252)];
tic
X1=dft(x,256);
toc
tic
X2=fft(x,256);
toc
时间已过 0.012424 秒。
时间已过 0.002118 秒。
x=[1,2,3,4,ones(1,1020)];
tic
X1=dft(x,1024);
toc
tic
X2=fft(x,1024);
toc
时间已过 0.249468 秒。
时间已过 0.008443 秒。
x=[1,2,3,4,ones(1,2044)];
tic
X1=dft(x,2048);
toc
tic
X2=fft(x,2048);
toc
时间已过 1.078018 秒。
时间已过 0.021545 秒
结果：
点数 DFT FFT
64 0.003188 秒 0.001967 秒
256 0.012424 秒 0.002118 秒
1024 0.249468 秒 0.008443 秒
2048 1.078018 秒 0.021545 秒
c
结果：
| 点数 | DFT | FFT |
|---|---|---|
| 64 | 0.003188 秒 | 0.001967 秒 |
| 256 | 0.012424 秒 | 0.002118 秒 |
| 1024 | 0.249468 秒 | 0.008443 秒 |
| 2048 | 1.078018 秒 | 0.021545 秒 |
分析：
5.利用FFT求信号频谱及分析采样频率、噪声对频谱的影响
考虑一个模拟信号源x(t)=3 \cos(8\pi t) + 6 \cos(20\pi t)按照时间间隔t=0.01n（其中n从0到N-1）进行采样和离散化处理，并绘制所得数字序列x_N(n)在不同点数（50和256点）下的DFT幅度谱图。通过快速傅里叶变换算法（FFT）进行计算和分析，并探讨频谱中所反映的数字频率域与实际频率之间的转换关系及其意义。
A程序流程如下：
- 初始化参数：设置采样周期为T_s = 0.01s
- 生成时间序列：根据给定公式生成连续时间函数值
- 进行数值积分：计算各时刻对应的离散值
- 应用DFT算法：利用快速傅里叶变换算法计算频域响应
- 绘制结果曲线：分别绘制50点和256点情况下的幅度分布
- 分析对比图形：比较不同点数对频谱显示效果的影响
N=50;
n=0:N-1;
T=0.01;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,1);
stem(n,x,'.');
title('f_s=100hz,N=50时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,2);
stem(n,abs(Y),'.');
title('f_s=100hz,N=50时信号的频谱');
N=256;
n=0:N-1;
T=0.01;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,3);
stem(n,x,'.');
title('f_s=100hz,N=256时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,4);
stem(n,abs(Y),'.');
title('f_s=100hz,N=256时信号的频谱')
c
结果：
B、采样频率为50hz，分别画出N=50，256点DFT幅度谱如下所示：
程序:
N=50;
n=0:N-1;
T=0.02;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,1);
stem(n,x,'.');
title('f_s=50hz,N=50时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,2);
stem(n,abs(Y),'.');
title('f_s=50hz,N=50时信号的频谱');
N=256
n=0:N-1;
T=0.02;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,3);
stem(n,x,'.');
title('f_s=50hz,N=256时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,4);
stem(n,abs(Y),'.');
title('f_s=50hz,N=256时信号的频谱');
c
结果：
分析：
程序：
将序列x_N(n)（其中N可选50或256，并需表示为一行向量）引入噪声0.2\cdot\text{randn}(1,N)；比较带与不带噪声时的信号频谱；增加强度至2\cdot\text{randn}(1,N)与$10\cdot\text{randn}(1,N)}；计算信噪比值；绘制并对比不同噪声水平下的时域波形图及幅频响应图；讨论不同强度下的结果对于信号分析的影响
figure(1)
N=50;
n=0:N-1;
t=0.01*n;
q=n*2/N;
x1=3*cos(8*pi*t)+6*cos(20*pi*t);
subplot(2,2,1)
plot(n,x1);title('Pure signal')
y=fft(x1,N);
subplot(2,2,2)
plot(q,abs(y));title('FFT of pure signal(N=50)')
x2=x1+0.2*randn(1,N);
SNR1=snr(x1,x2)
subplot(2,2,3)
plot(n,x2);title('noisy signal')
y=fft(x2,N);
subplot(2,2,4)
plot(q,abs(y));title('FFT of signal with noise (N=50)')
figure(2)
N=50;
n=0:N-1;
t=0.01*n;
q=n*2/N;
x1=3*cos(8*pi*t)+6*cos(20*pi*t);
x2=x1+2*randn(1,N);
x3=x1+10*randn(1,N);
SNR2=snr(x1,x2)
SNR3=snr(x1,x3)
subplot(2,2,1)
plot(n,x2);title('noisy signal（2*randn(1,N)）')
y2=fft(x2,N);
subplot(2,2,2)
plot(q,abs(y2));title('FFT of signal with more noise (N=50)')
subplot(2,2,3)
plot(n,x3);title('noisy signal（10*randn(1,N)）')
y3=fft(x3,N);
subplot(2,2,4)
plot(q,abs(y3));title('FFT of signal with more
noise(N=50)')
SNR1 =
0.0783
SNR2 =
-0.6612
SNR3 =
-6.4308
c
图形结果：
分析：
6.创新训练拓展内容
6.1 信号持续时间、频谱分析范围、采样点数和谱分辨率的关系。
已知模拟信号x(t)=\cos( 200\pi t)+\cos( 800\pi t)+\cos( 816\pi t)+\cos( 1000\pi t), 欲对其展开数字频谱分析, 并需确保能够清晰辨识出该信号的所有频率成分（任意两个相邻频谱分量之间至少间隔一个离散频率点）。为了便于采用基2-FFT算法进行频谱分析并避免频谱失真, 要求所选取的时间窗内必须包含该信号完整的整周期，并选取满足N=2M关系的采样点数（其中M为正整数）。请确定最低满足条件的采样频率（以避免发生频谱混叠）及其对应的最小记录时间, 并基于此设计并实现相应的数字系统硬件或软件平台以完成上述目标；最后完成相应的仿真实验并绘制出原始模拟信号的幅度频谱图, 分析所获得的频谱分析范围及所具有的分辨率。
通过计算可得f_{smin}=1024Hz, 对应所需的最小采样点数N_{min}=256个（具体推导过程见下文）。
基于上述参数设计相应的实验系统架构并完成实验验证工作流程。
T=1/1024;
n=0:255;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,1,1);
stem(n,x,'.')
title('采样1024hz的信号');
axis([0,256,0,4])
subplot(2,1,2);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=256;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*512,abs(X2))
xlabel('hz');
title('采样1024hz的幅度谱（N=256）');
axis([0,512,0,150])
c
结果：
分析：
（2）维持采样点数量不变，在采样频率上增加至原来的两倍后重新审视信号的幅度频谱特性。通过对比信号幅度频谱、以及相应的系统参数变化来探讨影响因素及其内在机理。
T=1/2048;
n=0:255;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,1,1);
stem(n,x,'.')
title('采样2048hz的信号');
axis([0,256,0,4])
subplot(2,1,2);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=256;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*1024,abs(X2))
xlabel('hz');
title('采样2048hz的幅度谱（N=256）');
axis([0,1024,0,180])
c
结果：
分析：
（3）维持采样频率恒定，在增加和减少一倍的采样点数下对信号进行重新分析。对比调整前后所得信号幅度频谱特征（包括谱分辨率和频谱分析范围），探讨这些变化所对应的系统响应特性差异及其影响因素
T=1/1024;
n=0:127;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,2,1);
stem(n,x,'.')
title('采样1024hz的信号');
axis([0,128,0,4])
subplot(2,2,2);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=128;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*512,abs(X2))
xlabel('hz');
title('采样1024hz的幅度谱（N=128）');
axis([0,512,0,85])
T=1/1024;
n=0:511;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,2,3);
stem(n,x,'.')
title('采样1024hz的信号');
axis([0,512,0,4])
subplot(2,2,4);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=512;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*512,abs(X2))
xlabel('hz');
title('采样1024hz的幅度谱（N=512）');
axis([0,512,0,300])
c
结果：
分析：
6.2频谱的内插函数恢复（频域采样恢复）
基于频域采样方法X(k)以及内插公式（教材式(3.3.7)和式(3.3.8)），通过Matlab编程实现离散频谱（DFT）向连续频谱的转换：在程序设计中选择N=4的矩形序列x(n)=R_4(n)，分别计算其4点和8点DFT值。采用内插公式恢复时，在Matlab程序中需确保ω的有效离散化。具体而言，在程序设计中要求至少使用200个离散点以确保足够的平滑度。随后生成X(e^{jω})的幅度频谱图，并对比不同DFT点数恢复效果的同时展示原始信号R_4(n)理论频谱。（参考教材图2.2.1）
解：程序：
subplot(2,2,1)
x=[1,1,1,1];
XK1=fft(x,4);
k=0:3;N=4;
w=[-2000:2000]*pi/1000;
X1=XK1/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi,abs(X1))
title('N=4的DFT内插公式恢复结果');
xlabel('w/π');
subplot(2,2,2)
n=0:3;
x=(1).^n;
w=[-500:500]*pi/250;
X=x*exp(-j*n'*w);
plot(w/pi,abs(X))
title('N=4的DFT的理论频谱');
xlabel('w/π');
subplot(2,2,3)
x=[1,1,1,1];
XK2=fft(x,8);
k=0:7;N=8;
w=[-2000:2000]*pi/1000;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi,abs(X2))
title('N=8的DFT内插公式恢复结果');
xlabel('w/π');
subplot(2,2,4)
n=0:3;
x=(1).^n
w=[-500:500]*pi/250;
X=x*exp(-j*(0:7)'*w);
plot(w/pi,abs(X))
title('N=8的DFT的理论频谱');
xlabel('w/π');
c
结果：
分析：
6.3对语音信号进行简单分析。
基于 MATLAB 平台设计开发了一套音频分析系统，在适当选取的采样频率下对不同个体采集到的单一音素语音信号进行预处理，并结合离散傅里叶变换方法实现对其频谱特性的研究工作。具体而言，在预处理时需对采集到的语音信号进行预处理以确保数据质量随后通过生成相应的时域波形图和频域 spectrogram 图直观展示其时序特性和频域特征进一步明确其对应的模拟频率范围并探讨其在不同频率范围内的分布特征对比不同个体在同一音素上的发音特征差异（建议选择具有代表性的男声与女声作为研究对象）并展开深入的技术探讨
figure(1)
recObj = audiorecorder;
disp('Start speaking.')
recordblocking(recObj, 2);
disp('End of Recording.');
play(recObj);
y = getaudiodata(recObj);
subplot(1,2,1);
stem(y,'.');
xlabel('时间');
ylabel('幅度');
title('信号波形');
subplot(1,2,2);
stem(abs(fft(y)),'.');
title('信号频谱');
xlabel('频率');
ylabel('幅度');
figure(2)
sec1=5
sec2 =12
y1=y(((1000*sec1+1):1000*sec2),:);
subplot(1,2,1);
stem(y1,'.');
xlabel('时间');
ylabel('幅度');
title('适当截取后信号的波形');
subplot(1,2,2);
stem(abs(fft(y1)),'.');
title('适当截取后信号的频谱');
xlabel('频率');
ylabel('幅度');
c
结果：
分析：
六、实验结论与心得体会（手写）
七、实验参考资料
高西全与丁玉美的著作《数字信号处理》由西安电子科技大学出版社在西安出版于2008年。
张德丰的《详解MATLAB数字信号处理》由电子工业出版社于北京出版于2010年。
王月明与张宝华合著的教材《MATLAB基础与应用教程》由北京大学出版社在北京出版于2012年。
陈怀琛编著的《数字信号处理教程：基于MATLAB的释义与实现》（第三版）由电子工业出版社在北京出版于2013年。
宋知用编著的《MATLAB数字信号处理85个实用案例精讲：从入门到进阶》由北京航空航天大学出版社在北京市出版于2016年。