电子电路]基于matlab的频谱分析及信号去噪仿真研究

电子电路]基于matlab的频谱分析及信号去噪仿真研究内容摘要：

)))*225。 figure。 imshow(A)。 在 MATLAB 中执行了 FFT 后，使用了 fftshift 函数调整，以使频谱图像的原点从起始点 (0， 0)，移到图像的中心点 (M/2,N/2)，对应的图 的傅立叶频谱能量图。 从图中可以看出：竖直方向理想单一频率的正弦分布纹理的频谱能量 集中在水平方向的三个点。 图 纹理频谱图 (3)离散信号 /序列 离散时间周期信号能够用具有谐波关系的复指数序列的线性组合来表示，称为离散傅里叶级数。 将这一概念推广应用到离散时间非周期信号，认为离散时间非周期信号也能够用具有谐波关系的复指数序列的线性组合来表示。 当离散时间周期信号的周期 N 趋于无穷大时，则离散时间周期信号就转化为离散时间非周期信号，其离散频谱就转化为连续频谱，称为离散时间傅里叶变换(Discrete Time Fourier Transform， DTFT)。 13 以 ( ) 1 2 sin ( 2 * 1 0 ) 3 sin ( 2 * 1 5 )x t t t  ，采样频率 100Hz，采样 100 个点，形成离散信号，利用 MATLAB 进行频谱分析。 程序如下： N=100。 Fs=100。 t=(0:(N1))/Fs。 xn=1+2*sin(2*pi*10*t)+3*sin(2*pi*15*t)。 XK=abs(fft(xn,N))。 f=(0:N/2)*Fs/N。 XK(1)=XK(1)/N。 XK(2:(N/2+1))=XK(2:(N/2+1))*2/N。 stem(f,XK(1:(N/2+1)))。 axis([1 N/2 0 5])。 grid on。 xlabel(39。 f(Hz)39。 )。 ylabel(39。 |X( f )|39。 )。 程序运行结果如图 34所示。 图 N=100 离散信号频谱 从图 可以看出， DFT 法分析的结果和实际信号吻合得很好，说明该方法确实有效。 但是，需要注意的是，离散傅里叶变换在频域是离散的，即限制在基频整数倍上，只能得到信号 K 次谐波成分的幅值谱，而对于非 K 次谐波成分的频谱则无法检测出来，并且由于栅栏效应和频谱泄漏，对其它 K次谐波的幅值也有影响。 如将上述程序中的 N=100 改为 N=50，结果如图 所示。 14 图 N=50离散信号频谱 可以看出，由于基频频率为 100/50=2Hz， x(t)中的 15Hz 分量未检测出来，其它频率分量的幅值也出现了误差，这种栅栏效应也是 DFT 应用中不可避免的问题之一。 要减小栅栏效应和频谱泄漏，可以减小采样频率，增加采样点数，以减小基频值，使谱线变密，这样原来漏掉的某些频谱成分就可以检测出来，但注意采样频率必须满足采样定理的要求，增加采样点数也会增加系统计算 DFT 的开销。 本章小结 在本章中， 前五节 介绍了频谱分析所用到的数学理论，分别为 采样定理，快速傅立叶变换，离散傅立叶变换， 介绍了频谱分析 基本原理以及 对数字信号处理的重要性， 阐述了频谱分析的方法， 在最后的 种 对几种 典型 信号进行了简单的频谱分析。 15 第 3 章 典型 信号去噪方法 IIR 和 FIR 滤波器的设计 IIR 滤波器和 FIR 滤波器的设计方法完全不同。 IIR 滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。 其设计步骤是：先设计过渡模拟滤波器得到系统函数 H（ s），然后将 H（ s）按某种方法转换成数字滤波器的系统函数 H(z)。 FIR 滤波器比鞥采用间接法，常用的方法有窗函数法、频率采样发和切比雪夫等波纹逼近法 [14]。 对于线性相位滤波器，经常采用 FIR 滤波器。 对于数字高通、带通滤波器的设计，通用方法为双线性变换法。 可以借助于模拟 滤 波器的频率转 换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。 具体设计步骤如下： （ 1）确定所需类型数字滤波器的技术指标。 （ 2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω =2/T tan( ) （ 3）将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。 （ 4）设计模拟低通滤波器。 （ 5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。 （ 6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。 我们知道 ，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。 为了克服之一缺点，可以采用双线性变换法。 下面我们总结一下利用模拟滤波器设计 IIR数字低通滤波器的步骤 [15] ： （ 1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频 率、阻带最小衰减。 （ 2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。 16 （ 3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。 （ 4）用双线性变换法， 模拟滤波器系统函数转换成数字低通滤波器系统函数。 如前所述， IIR 滤波器和 FIR 滤波器的设计方法有很大的区别。 下面我们着重介绍用窗函数法设计 FIR 滤波器的步骤。 如下： （ 1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角 窗 、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度 N。 先按照阻带衰减选择窗函数类型。 原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。 （ 2）构造希望逼近的频率响应函数。 （ 3）计算 h(n)。 （ 4）加窗得到设计结果。 在 Matlab 中，可以利用函数 fir1 设 计 FIR 滤波器，利用函数 butter,cheby1和 ellip 设计 IIR 滤波器，利用 Matlab 中的函数 freqz 画出各步步器的频率响应。 hn=fir1(M， wc， window)，可以指定窗函数向量 window。 如果缺省 window参数，则 FIR 默认为哈明窗。 其中可选 的 Rectangular Barlrtt Hamming Hann Blackman 窗，其相应的都有实现函数。 MATLAB 信号处理工具箱函数 buttp buttor butter 是巴特沃斯滤波器设计函数，其有 5 种调用格式， 本文中 用到的 是[N,wc]=butter(N,wc,Rp,As,’ s’ ),该格式用于计算巴特沃斯模拟滤波器的阶数 N和 3dB 截止频率 wc。 MATLAB 信号处理工具箱函数 cheblap,cheblord 和 cheeby1 是切比雪夫 I 型滤波器设计函数。 我们用到的是 cheeby1 函数，其调用格式如下：[B,A]=cheby1(N,Rp,wpo,’ ftypr’ ) [B,A]=cheby1(N,Rp,wpo,’ ftypr’ ,’ s’ ) 函数 butter,cheby1 和 ellip 设计 IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时 只需要代入相应的实现函数即可。 17 典型 信号的去噪对比 IIR 低通滤波器去噪 对正弦信号加噪声之后 再 进行去噪，利用 IIR滤波的方法 去除噪声。 正弦信号为 x=A*cos(2*pi*f*t/fs)， 其中： A=1， f=1000， fs=10^5。 程序代码如下： clf。 Ft=9000。 Fp=100。 Fs=300。 wp=2*pi*Fp/Ft。 ws=2*pi*Fs/Ft。 fp=2*Ft*tan(wp/2)。 fs=2*Fs*tan(wp/2)。 [n11,wn11]=buttord(wp,ws,1,50,39。 s39。 )。 %求低通滤波器的阶数和截止频率 [b11,a11]=butter(n11,wn11,39。 s39。 )。 %求 S域的频率响应的参数 [num11,den11]=bilinear(b11,a11,)。 %双线性变换实现 S域到 Z域的变换 [h,w]=freqz(num11,den11)。 %根据参数求出频率响应 plot(w*8000*,abs(h))。 legend(39。 用 butter设计 39。 )。 grid。 A=1。 %信号的幅值 f=1000。 %信号的频率 fs=10^5。 %采样频率 t=(0:999)。 %采样点 18 Mlag=100。 %相关函数长度变量 x=A*cos(2*pi*f*t/fs)。 %输入正弦波信号 xmean=mean(x)。 %正弦波信号均值 xvar=var(x,1)。 %正弦波信号方差 xn=awgn(x,5)。 %给正弦波信号加入信噪比为 20dB的高斯白噪声 figure(1) plot(t,xn)%绘制输入信号图像 title(39。 输入信号图像 39。 ) xlabel(39。 x轴单位： t/s39。 ,39。 color39。 ,39。 b39。 ) ylabel(39。 y轴单位： f/HZ39。 ,39。 color39。 ,39。 b39。 ) S=fft(x)。 %傅里叶变换 z11=filter(num11,den11,xn)。 m11=fft(z11)。 %求滤波后的信号 figure(1)。 subplot(2,1,1)。 plot(abs(S),39。 g39。 )。 title(39。 滤波前信号的频谱 39。 )。 grid。 figure(1)。 subplot(2,1,2)。 plot(abs(m11),39。 r39。 )。 title(39。 滤波后信号的频谱 39。 )。 grid。 figure(2)。 19 subplot(2,1,1)。 plot(xn)。 title(39。 滤波前信号的波形 39。 )。 grid。 figure(2)。 subplot(2,1,2)。 plot(z11)。 title(39。 滤波后的信号波形 39。 )。 0 100 200 300 400 500 600 700 800 900 10000200400600滤波前信号的频谱0 100 200 300 400 500 600 700 800 900 10000200400600滤波后信号的频谱 图 正弦 信号频谱图 图 为滤波 前 后的频谱对比图。 20 0 100 200 300 400 500 600 700 800 900 100042024滤波前信号的波形0 100 200 300 400 500 600 700 800 900 100021012滤波后的信号波形 图 去噪前后噪声对比 图 图 为滤波前后的信号 波形 对比图，可以看到，利用 IIR 低通滤波器虽然仍有小部分的噪声，信号并不完全平滑，但是绝大部分的 噪声已经被滤除。 FIR 低通滤波器去噪 设计 FIR 滤波器对信号进行去噪，信号为加噪声的正弦信号，信号表达式为x=cos(2*pi*f*t/fs)，其中， f=1000， fs=10^5。 程序代码如下： A=1。 %信号的幅值 f=1000。 %信号的频率 fs=10^5。 %采样频率 t=(0:999)。 %采样点 Mlag=100。 %相关函数长度变量 x=A*cos(2*pi*f*t/fs)。 %输入正弦波信号 xmean=mean(x)。 %正弦波信号均值 21 xvar=var(x,1)。 %正弦波信号方差 xn=awgn(x,5)。 %给正弦波信号加入信噪比为 20dB的高斯白噪声 figure(1) plot(t,xn)%绘制输入信号图像 title(39。 输入信号图像 39。 ) xlabel(39。 x轴单位： t/s39。 ,39。 color39。 ,39。 b39。 ) ylabel(39。 y轴单位： f/HZ39。 ,39。 color39。 ,39。 b39。 ) %fir滤波 wp=*pi。 %通带截止频率 ws=*pi。 %阻带截止频率 DB=wswp。 %过渡带宽度 N0=ceil(*pi/DB)。 M=N0+mod(N0+1,2)。 %计算 fir滤波器阶数 wc=(wp+ws)/2/pi。 %计算理想低通滤波器通带截止频率（关于π归一化） hn=fir1(M,wc)。 %调用 fir1计算 FIRDF的 h(n) y1n=filter(hn,1,xn)。 %将输入信号通过 fir滤波器 figur。