编程实现任意确定信号的频谱分析算法.doc

广西科技大学数字信号处理课程设计说明书题 目：编程实现任意确定信号的 频谱分析算法系 别：计算机工程学院专业班级：通信学 号：学生姓名：指导教师：目录摘要3一、设计内容3二、设计原理4三、设计过程71.CEG和弦音音频文件的频谱分析72.对该信号频谱能量较集中的频带滤波93.比较滤波前后的音频文件214.对滤波后的音频信号再滤出三个能量最集中的频簇215.重建信号与原信号的音频进行声音回放比较336.分析什么是和弦音35收获36参考文献36摘要：随着计算机和信息科学的飞速发展，信号处理逐渐发展成为一门独立的学科，成为信息科学的重要组成部分，在语音处理、雷达、图像处理、通信、生物医学工程等众多领域中得到广泛应应用。 Matlab语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，Matlab功能强大、简单易学、编程效率高，深受广大科技工作者的喜爱。特别是Matlab还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行语音信号分析、处理和设计。数字信号处理课程在现代科学中具有很大重要性及自身特点，理解与掌握课程中的基本概念、基本原理、基本分析方法，对用Matlab进行数字信号处理课程设计的思路，具有很大帮助。 语音信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。 关键词：抽样频率；频谱分析；滤波器 ；窗函数一、设计内容 （1）对给定的CEG和弦音音频文件取合适长度的采样记录点，然后进行频谱分析（信号的时域及幅频特性曲线要画出）。（2）分析CEG和弦音频谱特点，对该信号频谱能量相对较为集中的频带（分低、中、高频）实现滤波（分别使用低通，带通及高通），显示滤波后信号的时域和频域曲线，并对滤波后的信号与原信号的音频进行声音回放比较。（3）在低、中、高三个频带中，各滤出三个能量最集中的频簇，显示滤波后信号的时域和频域曲线。（4）任意选择几个滤出的频带（或频簇）进行时域信号重建（合成），与原信号的音频进行声音回放比较。讨论：根据上述结果，分析什么是和弦音二、设计原理数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。IIR滤波器的特征是，具有无限持续时间冲激响应。种滤波器一般需要用递归模型来实现，因而有时也称之为递归滤波器。FIR滤波器的冲激响应只能延续一定时间，在工程实际中可以采用递归的方式实现，也可以采用非递归的方式实现。数字滤波器的设计方法有多种，如双线性变换法、窗函数设计法、插值逼近法和Chebyshev逼近法等等。随着MATLAB软件尤其是MATLAB的信号处理工作箱的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化.线性相位FIR滤波器通常采用窗函数法设计。窗函数法设计FIR滤波器的基本思想是：根据给定的滤波器技术指标，选择滤波器长度N和窗函数（n），使其具有最窄宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h(n)。工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯塞（Kaiser）。FIR数字滤波器的设计原理：采用窗口法线性相位实系数FIR滤波器按其N值奇偶和h(n)的奇偶对称性分为四种：1、h(n)为偶对称，N为奇数 H(ej)的幅值关于=0，2成偶对称。 2、h(n)为偶对称，N为偶数 H(ej)的幅值关于=成奇对称，不适合作高通。3、h(n)为奇对称，N为奇数 H(ej)的幅值关于=0，2成奇对称，不适合作高通和低通。4、h(n)为奇对称，N为偶数 H(ej) =0、20，不适合作低通。 窗函数法设计线性相位FIR滤波器步骤: 确定数字滤波器的性能要求：临界频率k，滤波器单位脉冲响应长度N； 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ej)的幅频特性和相频特性； 求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ej)按M(M远大于N)点等距离采样，并对其求IDFT得hM(n)，用hM(n)代替hd(n)； 选择适当的窗函数w(n)，根据h(n)= hd(n)w(n)求所需设计的FIR滤波器单位脉冲响应； 求H(ej)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。 窗函数的傅式变换W(ej)的主瓣决定了H(ej)过渡带宽。W(ej)的旁瓣大小和多少决定了H(ej)在通带和阻带范围内波动幅度，常用的几种窗函数有：矩形窗 w(n)=RN(n)； Hanning窗 ； Hamming窗 ； Blackmen窗 ； Kaiser窗 。 式中Io(x)为零阶贝塞尔函数三、设计过程1.CEG和弦音音频文件的频谱分析对给定的CEG和弦音音频文件进行频谱分析，画出信号的时域及幅频特性曲线，程序如下：y,fs,bits=wavread(1.wav); %采集已录声音信号 N=length(y); %长度t=(0:length(y)-1)/fs; %计算时间figure;plot(t,y);xlabel(时间);ylabel(振幅); title(语音信号的时域图);Y=fft(y); %对信号进行FFTdf=fs/N; %采样间隔n1=0:N/2-1;f=n1*df; %频带宽度figure;plot(f,abs(Y(1:length(f);xlabel(频率);ylabel(振幅); title(语音信号的频谱图);sound(y); %播放给定的CEG和谐音频文件图形如下： 2.对该信号频谱能量较集中的频带滤波对该信号频谱能量相对较为集中的频带（分低、中、高频）实现滤波（分别使用低通，带通及高通）(1)FIR低通Ws=2*1200*pi/8000;Wp=2*1000*pi/8000;As=100;Rp=1; %数字低通滤波器参数Wc=(Wp+Ws)/2; %WcdeltaW=Ws-Wp; %过度带W取两个过渡带中的小者N0=ceil(12.8*pi/deltaW);N=N0+mod(N0+1,2); %根据表格7-3求出响应的窗口的点数N，保证N为奇数hd=ideallp(Wc,N); %低通滤波器的hd（n）%-ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录中，ideallp.mwdWindow=kaiser(N,10.056); %凯泽窗hr=hd.*wdWindow; %-点乘n=0:N-1; subplot(2,2,1); stem(n,wdWindow); %画二维离散图形 xlabel(时间);ylabel(振幅); title(凯泽窗); %x轴说明；y轴说明%标题subplot(2,2,2); stem(n,hr); %画二维离散图形xlabel(时间);ylabel(振幅); %x轴说明；y轴说明title(凯泽窗的冲激响应); %标题H,W=freqz(hr,1); %求每个窗对应的频率响应 subplot(2,2,3); plot(W/pi,abs(H); %画二维连续图形xlabel(omega/pi);ylabel(绝对振幅); %x轴说明；y轴说明title(滤波器低通幅频特性1); %标题legend(凯泽窗);subplot(2,2,4); plot(W/pi, 20*log10(abs(H); %画二维连续图形xlabel(omega/pi);ylabel(振幅(分贝); %x轴说明；y轴说明title(滤波器低通幅频特性2); %标题legend(凯泽窗);y,fs,bits=wavread(1.wav); %读出给定的CEG和谐音频信号N=length(y);t=(0:length(y)-1)/fs;df=fs/N;%采样间隔n1=0:N/2-1;f=n1*df;%频带宽度ly1=fftfilt(hr,y);%信号送入滤波器滤波，ly1为输出figure; %打开新的画图窗口subplot(2,2,1); plot(t,y);xlabel(时间); ylabel(振幅); title(语音信号的时域图);Y1=fft(y); %滤波前的FFT变换subplot(2,2,3);plot(f,abs(Y1(1:length(f); %画滤波前频谱图xlabel(频率); ylabel(振幅); title(语音信号的频谱图);grid;subplot(2,2,2);plot(t,ly1);xlabel(时间); ylabel(振幅);title(滤波后语音信号的时域图);Ly1=fft(ly1); %滤波后的FFT变换subplot(2,2,4);plot(f,abs(Ly1(1:length(f); %画滤波后频谱图xlabel(频率); ylabel(振幅);title(滤波后语音信号的频谱图);grid;wavwrite(ly1,fs,低通.wav);%把滤波后的信号写入文件sound(ly1); %回放低通滤波后的CEG和谐音频文件（2）FIR带通W1s=2*1000*pi/8000;W1p=2*1200*pi/8000;As=100;Rp=1;W2s=2*3200*pi/8000;W2p=2*3000*pi/8000; %数字带通滤波器参数W1c=(W1p+W1s)/2;W2c=(W2s+W2p)/2; %W1c,W2cdeltaW=min(W1p-W1s),(W2s-W2p); %过度带W取两个过渡带中的小者N0=ceil(12.8*pi/deltaW);N=N0+mod(N0+1,2); %根据表格7-3求出响应的窗口的点数N，保证N为奇数hd=ideallp(W2c,N)-ideallp(W1c,N); %带通滤波器的hd（n）wdWindow=kaiser(N,10.056); %凯泽窗hr=hd.*wdWindow; %-点乘n=0:N-1; subplot(2,2,1); stem(n,wdWindow); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴说明；y轴说明title(凯泽窗); %标题subplot(2,2,2); stem(n,hr); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴说明；y轴说明title(凯泽窗的冲激响应); %标题H,W=freqz(hr,1); %求每个窗对应的频率响应subplot(2,2,3); plot(W/pi,abs(H); %画二维连续图形xlabel(omega/pi);ylabel(绝对振幅); %x轴，y轴说明title(滤波器带通幅频特性1); %标题legend(凯泽窗);subplot(2,2,4); plot(W/pi,20*log10(abs(H); %画二维连续图形 xlabel(omega/pi);ylabel(振幅(分贝); %x轴，y轴说明title(滤波器带通幅频特性2); %标题legend(凯泽窗); y,fs,bits=wavread(1.wav); %读出给定的CEG音频文件N=length(y);t=(0:length(y)-1)/fs;df=fs/N;%采样间隔n1=0:N/2-1;f=n1*df;%频带宽度ly2=fftfilt(hr,y);%信号送入滤波器滤波，ly2为输出figure; %打开新的画图窗口subplot(2,2,1); plot(t,y);xlabel(时间); ylabel(振幅);title(语音信号的时域图);Y2=fft(y); %滤波前的FFT变换subplot(2,2,3);plot(f,abs(Y2(1:length(f); %画滤波前频谱图xlabel(频率); ylabel(振幅);title(语音信号的频谱图);grid;subplot(2,2,2);plot(t,ly2);xlabel(时间);ylabel(振幅); title(滤波后语音信号的时域图);Ly2=fft(ly2); %滤波后的FFT变换subplot(2,2,4);plot(f,abs(Ly2(1:length(f); %画滤波后频谱图xlabel(频率); ylabel(振幅);title(滤波后语音信号的频谱图);grid;wavwrite(ly2,fs,带通.wav);%把滤波后的信号写入文件sound(ly2,fs);%回放带通滤波后的CEG和谐音频文件（3）FIR高通Ws=2*3000*pi/8000;Wp=2*3200*pi/8000;As=100;Rp=1; %数字高通滤波器参数Wc=(Wp+Ws)/2; %WcdeltaW=Wp-Ws; %过度带W取两个过渡带中的小者N0=ceil(12.8*pi/deltaW);N=N0+mod(N0+1,2); %根据表格7-3求出响应的窗口的点数N，保证N为奇数 hd=ideallp(pi,N)-ideallp(Wc,N); %高通滤波器的hd（n）wdWindow=kaiser(N,10.056); %凯泽窗hr=hd.*wdWindow; %-点乘n=0:N-1; subplot(2,2,1); stem(n,wdWindow); %画二维离散图形xlabel(时间);ylabel(振幅); %x轴,y轴说明title(凯泽窗); %标题subplot(2,2,2); stem(n,hr); %画二维离散图形xlabel(时间);ylabel(振幅); %x轴.y轴说明title(凯泽窗的冲激响应); %标题H,W=freqz(hr,1); %求每个窗对应的频率响应 subplot(2,2,3); plot(W/pi,abs(H); %画二维连续图形xlabel(omega/pi);ylabel(绝对振幅); %x轴,y轴说明title(滤波器高通幅频特性1); %标题legend(凯泽窗);subplot(2,2,4); plot(W/pi,20*log10(abs(H); %画二维连续图形 xlabel(omega/pi);ylabel(振幅(分贝); %x轴,y轴说明title(滤波器高通幅频特性2); %标题legend(凯泽窗);%在一个窗口同时画出滤波前后的波形及频谱y,fs,bits=wavread(1.wav); %读出信号N=length(y);t=(0:length(y)-1)/fs;df=fs/N;%采样间隔n1=0:N/2-1;f=n1*df;%频带宽度ly3=fftfilt(hr,y);%信号送入滤波器滤波，ly3为输出figure; %打开新的画图窗口subplot(2,2,1); plot(t,y);xlabel(时间);ylabel(振幅);title(语音信号的时域图);Y3=fft(y); %滤波前的FFT变换subplot(2,2,3);plot(f,abs(Y3(1:length(f); %画滤波前频谱图xlabel(频率); ylabel(振幅); title(语音信号的频谱图);grid;subplot(2,2,2);plot(t,ly3);xlabel(时间); ylabel(振幅);title(滤波后语音信号的时域图);Ly3=fft(ly3); %滤波后的FFT变换subplot(2,2,4);plot(f,abs(Ly3(1:length(f); %画滤波后频谱图xlabel(频率); ylabel(振幅);title(滤波后语音信号的频谱图);grid;wavwrite(ly3,fs,高通.wav);%把滤波后的信号写入文件sound(ly3,fs);%回放高通滤波后的CEG音频文件3.比较滤波前后的的音频文件通过低、带、高通滤波后的音频文件与原音频文件有如下简单差异：低通滤波后，CEG和谐音频文件的声音较原来变得低沉了。带通滤波后，CEG和谐音频文件的声音没原来的那么尖锐。高通滤波后，CEG和谐音频文件的声音变得尖锐清亮。4对滤波后的音频信号再滤出三个能量最集中的频簇实现低、带、高通滤波后，在三个频带中，再分别用带通滤波器滤出三个能量最集中的频簇，显示滤波后信号的时域和频域曲线。（1）低通频簇W1s=2*100*pi/8000;W1p=2*200*pi/8000;As=100;Rp=1;W2s=2*900*pi/8000;W2p=2*800*pi/8000; %数字带通滤波器参数W1c=(W1p+W1s)/2;W2c=(W2s+W2p)/2; %W1c,W2cdeltaW=min(W1p-W1s),(W2s-W2p); %过度带W取两个过渡带中的小者N0=ceil(12.8*pi/deltaW);N=N0+mod(N0+1,2); %根据表格7-3求出响应的窗口的点数N，保证N为奇数hd=ideallp(W2c,N)-ideallp(W1c,N); %带通滤波器的hd（n）wdWindow=kaiser(N,10.056); %凯泽窗hr=hd.*wdWindow; %-点乘n=0:N-1; subplot(2,2,1); stem(n,wdWindow); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴，y轴说明title(凯泽窗); %标题subplot(2,2,2); stem(n,hr); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴，y轴说明title(凯泽窗的冲激响应); %标题H,W=freqz(hr,1); %求每个窗对应的频率响应subplot(2,2,3); plot(W/pi,abs(H); %画二维连续图形xlabel(omega/pi);ylabel(绝对振幅); %x轴，y轴说明title(滤波器带通幅频特性1); %标题legend(凯泽窗);subplot(2,2,4); plot(W/pi,20*log10(abs(H); %画二维连续图形 xlabel(omega/pi);ylabel(振幅(分贝); %x轴，y轴说明title(滤波器带通幅频特性2); %标题legend(凯泽窗);y1,fs,bits=wavread(低通.wav); %读出低通滤波后的信号t=(0:length(y1)-1)/fs;df=fs/N;%采样间隔n1=0:N/2-1;f=n1*df;%频带宽度ly11=fftfilt(hr,y1);%信号送入滤波器滤波，ly11为输出figure; %打开新的画图窗口subplot(2,2,1); plot(t,y1);xlabel(时间); ylabel(振幅);title(语音信号的时域图);Y11=fft(y1); %滤波前的FFT变换subplot(2,2,3);plot(f,abs(Y11(1:length(f); %画滤波前频谱图xlabel(频率); ylabel(振幅);title(语音信号的频谱图);grid;subplot(2,2,2);plot(t,ly11);xlabel(时间); ylabel(振幅);title(滤波后语音信号的时域图);Ly11=fft(ly11); %滤波后的FFT变换subplot(2,2,4);plot(f,abs(Ly11(1:length(f); %画滤波后频谱图xlabel(频率); ylabel(振幅);title(滤波后语音信号的频谱图);grid;wavwrite(ly11,fs,低通1.wav);%把滤波后的信号写入文件sound(ly11,fs);%回放（2）带通频簇W1s=2*1100*pi/8000;W1p=2*1200*pi/8000;As=100;Rp=1;W2s=2*2100*pi/8000;W2p=2*2000*pi/8000; %数字带通滤波器参数W1c=(W1p+W1s)/2;W2c=(W2s+W2p)/2; %W1c,W2cdeltaW=min(W1p-W1s),(W2s-W2p); %过度带W取两个过渡带中的小者N0=ceil(12.8*pi/deltaW);N=N0+mod(N0+1,2); %根据表格7-3求出响应的窗口的点数N，保证N为奇数hd=ideallp(W2c,N)-ideallp(W1c,N); %带通滤波器的hd（n）wdWindow=kaiser(N,10.056); %凯泽窗%-点乘hr=hd.*wdWindow;n=0:N-1; subplot(2,2,1); stem(n,wdWindow); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴，y轴说明title(凯泽窗); %标题subplot(2,2,2); stem(n,hr); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴，y轴说明title(凯泽窗的冲激响应); %标题H,W=freqz(hr,1); %求每个窗对应的频率响应subplot(2,2,3); plot(W/pi,abs(H); %画二维连续图形xlabel(omega/pi);ylabel(绝对振幅); %x轴，y轴说明title(滤波器带通幅频特性1); %标题legend(凯泽窗);subplot(2,2,4); plot(W/pi,20*log10(abs(H); %画二维连续图形 xlabel(omega/pi);ylabel(振幅(分贝); %x轴，y轴说明title(滤波器带通幅频特性2); %标题legend(凯泽窗);%在一个窗口同时画出滤波前后的波形及频谱y2,fs,bits=wavread(带通.wav); %读出带通滤波后的信号N=length(y2);t=(0:length(y2)-1)/fs;df=fs/N;%采样间隔n1=0:N/2-1;f=n1*df;%频带宽度ly22=fftfilt(hr,y2);%信号送入滤波器滤波，ly22为输出figure; %打开新的画图窗口subplot(2,2,1); plot(t,y2);xlabel(时间); ylabel(振幅);title(语音信号的时域图);Y22=fft(y2); %滤波前的FFT变换subplot(2,2,3);plot(f,abs(Y22(1:length(f); %画滤波前频谱图xlabel(频率); ylabel(振幅);title(语音信号的频谱图);grid;subplot(2,2,2);plot(t,ly22);xlabel(时间); ylabel(振幅);title(滤波后语音信号的时域图);Ly22=fft(ly22); %滤波后的FFT变换subplot(2,2,4);plot(f,abs(Ly22 (1:length(f); %画滤波后频谱图xlabel(频率); ylabel(振幅);title(滤波后语音信号的频谱图);grid;wavwrite(ly22,fs,带通1.wav);%把滤波后的信号写入文件sound(ly22,fs); %回放 （3）高通频簇W1s=2*3400*pi/8000;W1p=2*3500*pi/8000;As=100;Rp=1;W2s=2*4100*pi/8000;W2p=2*4000*pi/8000; %数字带通滤波器参数W1c=(W1p+W1s)/2;W2c=(W2s+W2p)/2; %W1c,W2cdeltaW=min(W1p-W1s),(W2s-W2p); %过度带W取两个过渡带中的小者N0=ceil(12.8*pi/deltaW);N=N0+mod(N0+1,2); %根据表格7-3求出响应的窗口的点数N，保证N为奇数hd=ideallp(W2c,N)-ideallp(W1c,N); %带通滤波器的hd（n）wdWindow=kaiser(N,10.056); %凯泽窗hr=hd.*wdWindow; %-点乘n=0:N-1; subplot(2,2,1); stem(n,wdWindow); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴，y轴说明title(凯泽窗); %标题subplot(2,2,2); stem(n,hr); %画二维离散图形 xlabel(时间);ylabel(振幅); %x轴，y轴说明title(凯泽窗的冲激响应); %标题H,W=freqz(hr,1); %求每个窗对应的频率响应subplot(2,2,3); plot(W/pi,abs(H); %画二维连续图形xlabel(omega/pi);ylabel(绝对振幅); %x轴，y轴说明title(滤波器带通幅频特性1); %标题legend(凯泽窗);subplot(2,2,4); plot(W/pi,20*log10(abs(H); %画二维连续图形 xlabel(omega/pi);ylabel(振幅(分贝); %x轴，y轴说明title(滤波器带通幅频特性2); %标题legend(凯泽窗);%在一个窗口同时画出滤波前后的波形及频谱y3,fs,bits=wavread(高通.wav); %读出高通滤波后的信号N=length(y3);t=(0:length(y3)-1)/fs;df=fs/N;%采样间隔n1=0:N/2-1;f=n1*df;%频带宽度ly33=fftfilt(hr,y3);%信号送入滤波器滤波，ly33为输出figure; %打开新的画图窗口subplot(2,2,1); plot(t,y3);xlabel(时间); ylabel(振幅);title(语音信号的时域图);Y33=fft(y3); %滤波前的FFT变换subplot(2,2,3);plot(f,abs(Y33(1:length(f); %画滤波前频谱图xlabel(频率); ylabel(振幅);title(语音信号的频谱图);grid;subplot(2,2,2);plot(t,ly33);xlabel(时间); ylabel(振幅);title(