数字信号处理实验报告

1、长春理工大学电工电子实验教学中心学生实验报告2014 2015 学年第一 学期实验课程数字信号处理实验实验地点 东4教学楼444实验室学院电子信息工稈学院专业通信工稈学号120421101姓绪杨杰实验项目音频信号釆集、分析及处理实验时间第十四周星期四9、10节实验台号操作成绩报告成绩三实验目的1. 以音频信号为例，熟悉模拟信号数字处理过程，进一步理解数字信号处理概念。2. 掌握运用matlab实现对音频信号的时频分析方法；3. 初步掌握数字音频信号合成的方法。4. 掌握运用matlab设计iir和fir滤波系统的方法；5. 掌握运用matlab实现对加噪音频信号进行去噪滤波的方法。锻炼学生运用

2、所学知识独立分析问题解决问题的能力，培养学生创新能力。二、实验仪器计算机、matlab软件三、实验性质综合分析、设计性实验四、实验任务实验内容一：windows系统中的“ding”音频信号的采集、分析、合成1. 音频信号的采集编写matlab程序，采集windows系统中的“ ding”声，得到*.wav音频文件，而后实现音频信 号回放。2. 音频信号的频谱分析运用matlab软件实现对音频信号的时域分析和频域分析，并打印相应的图形，完成在实验报 告中。注意：此音频信号的频谱包含两条主要谱线，在进行频谱分析时，注意频谱的完整性，利用matlab实现对两条主要谱线的定位并计算谱线所对应的模拟频率

3、。3. 音频信号的分解和合成运用matlab软件实现音频信号的分解与合成，将音频信号的频谱中两部分频谱成分进行分解，分别绘制出分解后的两个信号的频谱图；然后将分解后的两个信号再合成为一个新的信号，将合成后的新信号的时域、频域图与原来的信号时域、频域图相比较，绘制岀对比效果图。4. 音频信号的回放运用matlab软件实现音频信号的回放，将合成后的新信号和原音频信号分别进行回放，对比两个信号的声音效果。5. 音频信号分段傅里叶分析（选作）分析对一般音频.wav信号进行一次性傅里叶分析时存在的主要问题，利用分段傅里叶变换对该音频信号重新分析并合成。对比一次傅里叶分析结果并进行总结。实验内容二：任意音

4、频信号的时域和频域分析及数字滤波器设计1. 音频信号的采集音频信号的采集可以通过windows自带的录音机也可以用专用的录制软件录制一段音频信号（尽量保证无噪音、干扰小），也可以直接复制一段音频信号（时间为1s）,但必须保证音频信号保存为.wav的文件。2. 音频信号的时域、频域分析运用matlab软件实现对音频信号的打开操作以及时域分析和频域分析，并画出相应的图形，打印在实验报告中。3. 引入干扰信号在原有的音频信号上，叠加一个频率为100khz的正弦波干扰信号（幅度自定，可根据音频信号幅度情况而定）。4. 数字滤波器设计运用matlab实现iir或fir滤波系统，要求加入干扰的音频信号经过

5、滤波系统后，能够滤除 100khz的干扰信号，同时保留原有的音频信号，要求绘制出滤波系统的频谱图，同时绘制出滤波 前后音频信号的时域、频域对比图。实验内容三：编写 fft算法程序（此部分为选作内容）分析fft算法实现过程，运用 matlab软件实现fft函数功能，要求运用自行编写的matlab程序实现对固定信号（自定义信号）的频谱分析，并绘制出频谱图；然后将自己设计的fft算法程序与fft函数实现的频谱图相比较，分析两个程序结果的差别及原因。五、实验任务分解及思路实验内容一：windows系统中的“ding”音频信号的釆集、分析、合成通过对任务内容要求的分析，建立程序框图如图一（任务一程序框图

6、）所示，详细程序及结果见附录。编程思想及思路如下：1、首先用audioread读取音频信号得到 signal序列及采样频率fs,由于音频信号是双声道，这里只取第一个声道进行处理，由于原信号本身就较短，这里对长度就不作处理。接着按实验要求进行音频回放，并绘制出音频信号的时域波形。2、调用fft函数，对signal进行快速傅里叶变换，并绘制出频谱图。在绘制频谱图时考虑到快 速傅里叶变换的对称性，这里只取序列的前半部分进行观察分析。3、为了实现音频信号的分解及合成，先对原信号的频谱图进行观察分析，发现原信号的主要能量集中在两个频率上，于是考虑用这两个频率的正弦信号合成该信号。这样首先要找到 这两个频

7、率。为了求得这两个频率，我先调用findpeaks函数找到频谱图上的各个峰值点signal_pks及其对应的位置iocs ,然后用sort对峰值点进行排序，找到最大的两个值， 接着用find函数找到这两个峰值点在iocs中的位置，这样就从iocs中知道了最大两个峰值点在原音频信号中的位置，也就知道了对应的频率。接着生成这两个频率的正弦信号，幅度的比例等于频谱图上对应频率的幅值比例。然后将两个正弦信号叠加，得到合成后的信号，绘制出时域波形，与原信号波形进行比较，接着对两个正弦信号进行fft ,绘制岀他们的频谱，然后对合成的信号进行fft ,做出频谱图与原信号的频谱图进行比较。4、用sound函数

8、进行原信号和合成信号的回放，比较差异。5、要对音频信号进行分段傅里叶分析，首先就要确定每段的长度。考虑计算机进行傅里叶变换时点数为时，为便于计算，这里就将每段长度暂定为1024,在for循环中对每段信号用fft函数进行快速傅里叶变换，然后再用ifft对每段快速傅里叶变换后的序列进行快速傅里叶逆变换并合成，然后进行音频回放，与原信号进行比较。寻找频谱最高护个峰值点找到峰值点苓应的頻率等于原信号最大幅值做合成信号时域波 些原罕出交 做合成信号丁域旅丈图 士故两旺弦罕的频谱做合成信号的频谱与原信吗華i普比孙原信号及合成信号的回护较原信号的分段傅里图一任务一程序框图实验内容二：任意音频信号的时域和频域

9、分析及数字滤波器设计通过对任务内容要求的分析，建立程序框图如图二（任务二程序框图）所示，详细程序及结果见附录。编程思想及思路如下：1、首先用audioread读取音频信号signal及采样频率fs,同实验内容样只取第一声道信号2的数据进行处理。进行分析处理，由于此处选择的音频文件较长，这里截取前2、同实验任务一一样，做出信号的时域波形，及 fft变换后的频谱图。图二任务二程序框图3、这里要求引入100khz的正弦干扰信号，由于由 1中已得到fs为22050h乙 根据奈奎斯特频率釆样定理，采样频率必须大于等于原信号最高频率的2倍，所以必须将原信号的采样频率提高到200khz以上才能引入100kh

10、z的噪声，所以这里考虑用一阶线型插值interpl将原信号的采样频率提高到220500h乙这样就可以引入 100khz噪声。做出提高采样频率后的信号的时域波形和频谱图，确认信号并没有发生变化。接着生成100khz的正弦信号，根据2中做出的信号的时域波形的幅度，这里取噪声的幅值为0.5o将提高采样频率后的4、信号与噪声叠加。对加噪后的信号做出时域波形和频谱图，观察波形的变化。这里要求设计数字滤波器，根据对加噪前的频谱以及加噪后的频谱的观察，可以采用低通考虑设计相对较为简单的巴特沃斯低通滤波器进行滤波。滤波前首先要确定设计指标，观察频谱这里暂取0然后开始设计巴特沃斯低通滤波器,这里我把设计的程序打

11、包成一个函数方便调用，函数的框图如图三（巴特沃斯低通滤波器开始k= 1 , n参量输入qc = ftp(10°lap-1)nc = np(io°lap 由rz,pz,kz合幷部分分式得hr(z)由rz,pz,kz合幷部分分式得hr(z)2avs(k) = nrp(fc)jf«(s项式合并确定九分子多项式民(s)部分分式展开得rs.psksz(k)=訂只如rss=rs/fsx;>pz(k.l)=iyz(k)nlength b < n+1| b(k) = 0(*!=；bz =bz= mzla| az=flip(a)n输出占乙/结束图三 巴特沃斯低通滤波器函

12、数运算流图)可以发现fft的运算过程很有规律性。h二2"点的m共进行m级运算，每级由n/2个a(0)a(0)a(0)a(0)x(0)x(1)xxxxxx8点dit-fft运算流图图四函数)所示。该函数调用需要5个输入参数分别为：n.：模拟通带边界角频率，qs：模拟 阻带边界角频率，耳：通带最大衰减，l：阻带最小衰减，fs：信号釆样频率。参数输入后，根据数学上的推导开始计算。首先根据书本公式6.2.18确定滤波器的阶数 n,然后根据书本公式6213求出归一化极点卩卜 接着求出3db截止频率qy然后对卩£进行去归一化 得到» 然后得到h/s)的分子分母系数，接着对进行

13、部分分式展开，由书本公式6.3.6, 利用脉冲响应不变法将f4(s)变换到h(z)，接着合并部分分式，求得 h(z)的分子分母多项式的系数，对多项式的系数进行整理排列得到巴特沃斯低通滤波器系统函数的分子分母多项式的系数向量b乙az,作为函数输出。接着调用freqz函数求出设计的低通滤波器的频率 响应特性函数，然后绘制出去频响特性曲线并观察，看是否符合设计指标。若符合则开始 滤波，不符合就要重新设计。滤波时直接调用了 filter函数，得到的信号做出时域波形和频 谱图并与滤波前的信号对比，观察滤波效果，是否滤除了100khz的干扰信号，并保留了原有的音频信号。实验内容三：编写 fft算法程序此部

14、分要求自己编写fft算法程序，首先要充分了解fft的运算过程。观察图四(8点dit-fft蝶形运算组成。同一级中每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元。这样，经过 m级运算后，原来存放输入序列数据的n个存储单元(数组 a)中便一次存放 x(k)的n个值。这样可以节省大量内存，从而降低设备成本。转因子的指数。但各级的旋转因子和循环方式都有所不同。为了编写计算程序，应先找出旋转因子与运算级数的关系。用表示从左到右的运算级数。 不难发现，第l级共有p :个不同的旋转因子。

15、对一般情读入x(n),m况，第 l级的旋转因子为昭遡斛冃仏严因为2、山尸二忆曲所这样就确定了第 l级运算的旋转因子。接着研究蝶形运算规律，设序列 x(n)经吋域抽选(倒序)后，按图中次序存入数组a中，如果蝶形运算的两个输入数据相距b个点，应用原位计算，则蝶形运算相距b个点，应用原位计算，则蝶形运算开始倒序j= o.b-1表示成如下形式心(/)"“(/)+也/+3)昭 al(j+b)t = a(k)+a(lc + b)昭/ = 0x2.m,2w-1；下标l表示第l级运算，4(/)则表示第l级运算后的数组元素 a(j)的值。ffu-10)表示第l输出结束级运算前a(j)的值。图五 fft

16、运算程序框图另外由图四还有一些运算规律：第l级中，每个蝶形的两个输入数据相距3二?个点；每级有b个不同的旋转因子；同一旋转因子对应着间隔为卩点的个蝶形。总结上述规律，采用以下算法。先从输入端开始，逐级进行，共进行m级运算。在进行第 l级运算时，依次求岀 b个不个蝶形。这样，利用三重循环实同的旋转因子，每求出一个旋转因子，就计算完它对应的所有 现fft算法，程序框图如图五(fft运算程序框图)所示。忙外，ditffti云算流图的输出 x(k)为白然顺丿命，但一为了i舌应原付计算，其输入丿命列不是按的自然顺序排列，憩建m次偶奇抽选后的排序称为序列x(n)的倒序。因此，在运算 m级蝶形之前应对序列x

17、(n)进行倒序。由于本程序使用的是 matlab高级语言无法直接倒置二进制数位， 因 此必须找到产生倒序的十进制运算规律。观察倒序的二进制数可以发现，倒序数是在m位二进制数最高位加1,逢2向低位进位。为了叙述方便，用j表示当前倒序数的十进制数值。对于h二2" m位二进制数最高位的十进制权值为n/2,lh=n/2j =lh2 =n-21且从左向右二进制位的权值依次为n/4, n/8,，2, 1o因此，最高位加1相当于十进制运算 j+n/2。如果最高位是0 (jvn/2),则直接由j+n/2得下一个倒序值；如果最高位是1 ( j> n/2),则先将最高位变成 。问7/2),i &g

18、t; j ?ni = 1 , n1然后次高位加 1(j+n/4) o但次高位加1时，同样要断0、1值；依次类推直到完成最高位加逢2向右进位的运算。倒序的程序框图如图六 (倒序程序框图)所示。形成倒序j后，将元数组中存放的输入序列重新按倒序排序。 设原输入序列x(l)先按自然顺序存放入数组n?j = j-kk = k/2t = a(l) a(l) = a(j)进行倒序，发现第一个勺最后一个虚列值不需要重排， 所以顺序数图六倒序程序框图中。对x(n)i的初值为1,盤为n-2,倒序数j的初值为n/2o每计算出一个倒序值j便与循环语句自动生成的顺序i比较,当i二j吋，不需要交换，当 冈吋，a和a(j)

19、交换数据。另外，为了避免重复调换，对ivj的情况进行调换。以上是fft函数的编程思路，具体程序见附录，已純yj_fft(函数，方便调用，附录中遞对同一个音频序列用 yj_fft(和fft()进行傅里叶变换的比较测试程序，用于观察设计的函数的数 性。八、实验任务windows系统中的“ding喑频信号的采集、分析、合成(程序框图见图一）contents.音频信号的采集12音频信号的频谱分析自动搜寻频谱的峰值点并在图中表示出来1.3.音频信号的分解和合成1.4音频信号的回放1.5音频信号分段傅里叶分析% % title:数字信号处理实验一一音频信号采集、分析及处理% author : 120421

20、101 杨杰% % % 实验任务一：windows系统中的“ding”音频信号的采集、分析、合成% %clear all ;close all ;1.1 音频信号的采集sign al,fs = audioread('ding.wav* );淤取音频文件，得到音频序列signal及采样频率fssignal = sign al(:,1);%只取第一个声道的序列len_signal = length signal);% 言号的长度t= 1/fs;嫁样周期t = t:t:le n_ signal/fs;%!寸间t序列sound( signal,fs);%原音频信号播放pause(1);痼停1s

21、1.2.音频信号的频谱分析figure 'name*音频信号的频谱分析'position ,50,350,400,300）;涯成一个窗口subplot( 2,1,1 );刘各窗口分为两行，第一部分用于绘制时域波形plot( t,signal);滋制时域波形axis( t len_signal/fs -0.08 0.08);%艮定显示范围title（|时域波形图'）;淞制该波形的标题xlabel( 1 时间(s); ylabel(1 幅值')；滋制横纵坐标的标签fft_n = 2anextpow2( len_signal);加十算fft变换的点数signal =

22、fft( signal,ffi_n );%快速傅里叶变换signal.amp = abs( signal );%只做幅频分析，对fft变换后序列取模f = fs*(0:fft_n/2-1)/fft_n;瀰率横轴f,由于对称性，只取前一半signal amp h/lfv4gn/l amp (斑%subplot( 2,1,2 );%窗口第二部分，绘制频谱图08唸会制频谱图唏示题囁 会制横纵坐标的标签哪艮定显示范围plot( f,signal_amp_half );title('频谱图)；xlabel( * 频率(hz );ylabel( * 幅度); axis( 0 4000 0 250)

23、;时域波形图0 050-0.050406时间 频谱图b频率(hz)下面程序对信号的分解分析：观察音频信号的频谱图可知该音频信号能量主要集中在两个频率上, 主要是找到这两个峰值点对应的频率。自动搜寻频谱的峰值点并在图中表示出来signal.pksjocs = findpeaks( signal_amp_half ); sort_pks = sort( signal_pks, 'descend'); sortjoc = find( signal_pks>=sort_pks(2);rd 二 locs( sort_loc(1);n2 = locs( sort_loc(2);y1

24、= signal_amp_half(n1);y2 = signal_amp_half(n2);f1 = (n1-1)*fs/fft_n;f2 = (n2-1)*fs/fft_n;s1 = sprintf(* 峰值一 (%fhz)' ,f1 );s2 = sprintf(* 峰值二(%fhz)' ,f2);text( f1,y 1,s1 );text( f2,y2,s2 );%寻找频谱峰值点加寸频谱峰值点进行排序呀戈到峰值最大的2个点排序后在pks中的位置艰高峰值点在序列中的位置癖二高峰值点在序列中的位置艰高峰值点在频谱中的幅度瞬二高峰值点在频谱中的幅度鴻高峰值点频率獅二高峰值点

25、频率漩义字符串1漩义字符串2陆图上最高峰值点处显示“峰值一”拓图上第二高峰值点处显示“峰值二”时域波形图0.05 0-0.05020406时间0.8频率(hz)分析：找到了两个峰值点对应的频率并在图上显示出来。1.3.音频信号的分解和合成y1 =(y1/(y1+ y2)*cos(2*pi*f1*t); y2 = (y2/(y1+ y2)*cos(2*pi*f2*t);瞬一个频率信号瞬二个频率信号a = max( signal);hecheng = (y1 +y2)/max( y1+y2 );hecheng = hecheng*a;%原信号的最大幅度殆成信号归一化酬吏合成信号与原信号的幅度相等f

26、igure( name*音频信号的合成，position1 subplot( 2,1,2 );plot( t,hecheng );title( *合成音频时域波形图)；xlabel( * 时间(s) ); ylabel，幅值)；,450,350,400,300);傩成一个窗口加各窗口分为两行，第二部分用于绘制时域波形呀会制合成信号时域波形唏会制该波形的标题淤制横纵坐标的标签subplot( 2,1,1 );h1=plot( t,signal);axis( t len_signal*t -0.08 0.08); title('原信号时域波形图1 );xlabel( * 时间(s) );

27、ylabel( 1 幅值);煖口第一部分，绘制原信号时域波形喏会制原信号时域波形哪艮定显示范围滋制标题淞制横纵坐标的标签figure( 'name* ,'合成音频时域波形图(时域放大) plot( t,hecheng );title('合成音频时域波形图(时域放大)；xlabel('时间(s)'); ylabel( '幅值)；axis 0.2 0.21 -0.08 0.08);1 , position* ,850,350,400,300);成一个窗口淞制时域波形淞制标题淞制横纵坐标的标签滋艮定显示范围(局部放大)%生成一个窗口figure( &#

28、39;name'， 音频信号的分解,'position1 ,50,10,400,300);subplot( 2,1,1 );fft_y1 = abs( fft( y1*a,fft_n );plot( f,fft_y11:(fft_n/2);title('第一个频谱峰值点信号频率谱)；xlabel( * 频率(hz);text( f1,fft_y1(n1),s1 );axis 0 4000 0 1000);subplot( 2,1,2 );fft_y2 = abs( fft( y2*a,fft_n );plot( f,fft_y2(1:(fft_n/2);title(&#

29、39;第二个频谱峰值点信号频率谱)；xlabel( * 频率(hz );ylabel( * 幅度);text( f2,fft_y2(n2),s2);axis( 0 4000 0 1000);figure( name*合成前后频率谱比较，position* subplot( 2,1,2 );hecheng = abs( fft(hecheng,fft_n);plot( f,hecheng(1:(fft_n/2);title('合成后信号频率谱-)；xlabel( 1 频率(hz );ylabel'幅度)；text( f1,fft_y1(n1),s1 );text( f2,fft_

30、y2(n2),s2);axis 0 4000 0 1000);subplot( 2,1,1 );plot( f,signal_amp(1:(fft_n/2);title('原信号频谱1 );xlabel('频率(hz );ylabel'幅度)；axis 0 4000 0 250);text( f1,y1,s1 );text( f2,y2,s2 );剌各窗口分为两行，第一部分用于绘制第一个频率信号的频谱%对第一个频率信号进行 fft ,并取模 您制频谱图%绘制标题淹制横纵坐标的标签%在峰值点处显示“峰值一” 卅艮定显示范围爛口第二部分，第二个频率信号的频谱 慟第二个频率信

31、号进行fft ,并取模 淞制频谱图%会制标题滋制横纵坐标的标签独峰值点处显示“峰值二” 哪艮定显示范围,450,10,400,300);淞成一个窗口爛口第二部分，绘制合成信号频谱图 腋寸合成信号进行fft变换并取模 燈制频谱图淞制标题潍制横纵坐标的标签峰值点处显示“峰值一” 憾峰值点处显示“峰值二” 滋艮定显示范围嘀口第一部分，绘制原信号频谱图 淞制原信号频谱图滋制标题滋制横纵坐标的标签%艮定显示范围潍图上最高峰值点处显示“峰值一” 注图上第二高峰值点处显示“峰值二”原信号时域波形图时间合成音频时域波形图010.2040.60.8时间.20 2020 2040 2060 208时间合成音頻时域

32、波形图（时域放大）0.080.060.040.02凹 00.02-0.04-0.060峠分析：观察音频信号合成后的时域波形可以发现合成后的音频信号与原始音频信号差距较大，这是因为原始音频信号从频谱图上看，虽然看似能量都集中在两个频率上，但实际上，在那两个频率的 两旁还有许多频率分量，用这种音频信号合成方法合成音频信号时，忽略掉了其他的所有频率分量， 合成出来的音频信号自然与原信号有很大差异。第一个频谥峰值点信号频率谱顎率(hz)第二个频语峰值点信号频率谱1000 黒 500梅值二(1242 869568hz)-90 1 01000200030004000频率(hz)原信号频谱200磬100峰値

33、一(93131035hz) 殍值二(1242 869568hz)°01000100020003000频率(hz)合成后信号频率谱予伯35甩)_l4000值二(1242 869568hz)°0100020003000频率(hz)4000%原音频信号播放瞬停1s殆成音频信号播放瞬停1s可以发现分析这是分解的音频信号的频谱图以及合成后的频谱图和原信号的频谱图的比较。计算得到的两个频谱峰值对应的频率是正确的。1.4音频信号的回放sound( signal,fs); pause(1);sound( hecheng,fs); pause(1);1.5音频信号分段傅里叶分析n = 102

34、4;m = ceil( len_signal/n );fftfenduan 二 zeros( m,n );for k =fftfenduanf k,:) = fft( siqnal(n*(k-1)+1:n*k),n );滋定每段的点数2切，此处为1024漏定分的段数炬成n行m列数组用于存储分段fft后的序列%分段傅里叶变换endfftfenduan( m,:) = fft( signal(n*k+1:len_signal),n );fenduan = zeros( 1,m*n );for k = 1:m;fenduan(1,n*(k-1)+1:n*k) = ifft( fftfenduan(k

35、,:),n );end演后一行音频信号进行 m%生成m5l长序列，用于存储ifft 后的序列%分段傅里叶逆变换并合成%分段傅里叶分析后的音频信号的播放,850,10,400,300);炬成一个窗口頰口第一部分，绘制分段傅里叶分析后信号时域波形 %会制分段傅里叶分析后信号时域波形哪艮定显示范围強会制标题唸会制横纵坐标的标签滋寸分段傅里叶分析后信号做fft ,并取模煖口第二部分，绘制分段傅里叶分析后信号频谱图強会制分段傅里叶分析后信号频谱图唸会制标题您制横纵坐标的标签哪艮定显示范围sound( fenduan,fs );figure( name，,'分段傅里叶分析后时域频域图，positi

36、on*subplot(2,1,1);plot( t,fenduan(1:len_signal);axis t len_signal/fs -0.08 0.08);title(*分段傅里叶分析后音频时域波形图)；xlabel( * 时间(s); ylabel( * 幅值)；%分段分析合成后的傅里叶变换频谱fenduan = abs( fft( fenduan,fft_n );subplot(2,1,2);plot( f,fenduan(1:floor(fft_n/2);title('分段分析后音频频率谱* );xlabel('频率(hz );ylabel'幅度)；axis

37、( 0 4000 0 250);分段傅里叶分析后音频时域波形图0 050-0.050.40.6时间08频率(hz)分析：这是分段傅里叶分析后合成的音频信号的时域波形及频谱图，发现和原始信号一样，相对于 一次性傅里叶分析，很好地实现了信号的合成。因为一次性傅里叶分析过程中只考虑了单一的两个 频率分量，忽略了其他所有的频率分量，而分段傅里叶分析考虑了所有的频率。总结，对于音乐或 说话的声音，一般不可能只有单一的几个频率合成，它们往往是由许许多多的频率分量组成的，考 虑的频率越多，所合成的信号越接近原始的音频信号。published with matlab? r2014a实验内容二：任意音频信号的时

38、域和频域分析及数字滤波器设计（程序框 图见图二）contents2.1 音频信号的采集22音频信号的频谱分析23引入干扰噪声100k2.4.数字滤波器设计25用数字滤波器进行滤波调用flite函数% title :数字信号处理实验一一音频信号采集、分析及处理% author : 120421101 杨杰 % % % 实验任务二：任意音频信号的时域和频域分析及数字滤波器设计% %clear all ;close all ;2.1 音频信号的采集sign al,fs = audioread('reny i.wav1);%读取音频文件，得到音频序列signal及采样频率fslen_sig n

39、al = fs;signal = signal(1:le resign al,1);%截取信号的长度为1s%只取第一个声道的序列t= 1/fs;殊样周期t = t:t:len_s ign al*t;%时间序列t,50,350,400,300);2.2.音频信号的频谱分析figure( 'name'， 音频信号分析','position* subplot( 2,1,1 );plot( t,signal);title('时域波形图* );xlabel( 1 时间,);ylabel( 1 幅值')；fft_n = 2anextpow2( len_sig

40、nal);signal = fft( signal,fft_n );signal.amp = abs( signal );f = fs*(0:fft_n/2-1)/fft_n;signal_half 二 signal_amp(1:(fft_n/2);%生成一个窗口%各窗口分为两行，第一部分用于绘制时域波形%绘制时域波形%绘制该波形的标题%绘制横纵坐标的标签硏算fft变换的点数漱速傅里叶变换頫做幅频分析，对ffi变换后序列取模獅率横轴f,由于对称性，只取前一半%由于对称性丿貝取前一半subplot( 2,1,2 );滴口第二部分，绘制频谱图嗨会制频谱图術题滋制横纵坐标的标签%艮定显示范围%提高信

41、号的采样频率10倍，以满足可以引入fs = fs*10;tt= 1/fs;tt = t:tt:len_signal*t;signalfs = interpl (t,signal,tt); len_signalfs = length( signalfs );figure( name，,'提高采样频率后的音频信号 subplot( 2,1,1 );plot( tt,signalfs);title(1提高采样频率后音频时域波形图xlabel( * 时间(s); ylabel( * 幅值)；100khz的噪声%各采样频率提高10倍 淞高采样频率后的采样周期 加寸间序列tt瀏原信号一阶线型插值

42、漩值后的音频信号的长度 , position* ,450,350,400,300);池成一个窗口剌各窗口分为两行，第一部分用于绘制采样频率提高后的时域波形滋制釆样频率提高后的时域波形嚇题強会制横纵坐标的标签fft_nn = 2anextpow2( len_signalfs );signalfs = fft( signalfs,fft_nn );signalfs.amp = abs( signalfs );ff = fs*(0:fft_nn/2-1)/fft_nn;signalfs_amp_half = signalfs_amp( 1:(fft_nn/2);subplot( 2,1,2 );pl

43、ot( ff,signalfs_amp_half );title('提高采样频率后音频频谱)；xlabel( * 频率(hz );ylabel( * 幅度);axis( 0 5000 0 4000);力十算采样频率提高后fft变换的点数 於夬速傅里叶变换頫做幅频分析，对fft变换后序列取模 獅率横轴ff ,由于对称性，只取前一半 %由于对称性，只取前一半癖二部分用于绘制采样频率提高后的频谱图滋制采样频率提高后的频谱图嚇题淞制横纵坐标的标签加艮定显示范围0020.40.6081时间频谱plot( f,signal_half );title('频谱)；xlabel( * 频率(hz

44、 );ylabel( * 幅度); axis( 0 5000 0 400);频率(hz)0.2提高采样頻率后音频时域汲形團-0.20020.406081时间提高采祥频率后音频频谙劭率(hz)分析：首先看原信号的时域波形及频谱图，可知该信号截取的是 1s,频率分量主要集中在低频段, 提高采样频率后比较时域波形及频谱图还是和原始信号一样，说明提高采样频率过程中并没有改变 信号，只是提高了采样频率，符合要求。2.3.引入干扰噪声100knoise = 0.5*cos(2*pi*100*10a3*tt);signal_noise = sign alfs + no ise1:le ngth(signal

45、fs); figure( 'name* , 1引入噪声后的音频信号,'position* subplot( 2,1,1 );plot( tt,signal_noise);title('引入噪声后音频时域波形图)；xlabel( * 时间(s) ); ylabel，幅值)；%100khz 噪声啕入噪声,850,350,400,300);淞成一个窗口刘各窗口分为两行，第一部分用于绘制引入噪声后的时域波形 啕【入噪声后的时域波形 滋示题淤制横纵坐标的标签加寸加噪后信号进行fft頫做幅频分析，对fft 变换后序列取模%由于对称性，只取前一半瞬二部分用于绘制加噪声后音频频率谱 淞

46、制加噪声后音频频率谱 術题滋制横纵坐标的标签signal.noise = fft( signal_noise,fft_nn );signal_noise_amp = abs( signal.noise );signalfs_noise_amp_half = signal_noise_amp( 1:(fft_nn/2); subplot( 2,1,2 );plot( ff,signalfs_noise_amp_half );title('加噪声后音频频谱1 );xlabel( * 频率(hz)，);ylabel( * 幅度);引入哩声后音频时域波形圉liiii"1o021020

47、406时间加噪声后音频倾谱0.0°0468頻率(hz)10 12x1owp = 2*pi*4000;ws = 2*pi*8000;ap = 1；as = 60;b乙az = yj_butter_low(wp,ws,ap,as,fs);h,w = freqz( b乙az);hf = abs( h );figure( name，,'巴特沃斯低通滤波器系统函数图 subplot( 2,1,1 );plot( w,hf);title('滤波系统幅频特性曲线* );xlabel('数字频率(0兀);ylabel( * 幅度,); subplot( 2,1,2);plot

48、( w,20*log10(hf);title('滤波系统幅频特性曲线(db) * ); xlabel(数字频率(0龙)*); ylabel( 'db');滅离散系统频率响应特性函数 颁模分析：这是加噪后的信号时域波形及频谱图，可以知道该信号时域波形发生了较大的改变，从频谱 图中也可以发现在10人5处即1ookhz处出现了一个较高的冲激。说明成功引入了 1ookhz的正弦干 扰信号。2.4.数字滤波器设计%巴特沃斯低通滤波器带边界角频率滋且带边界角频率带最大衰减% 且带最大衰减澜自己设计的巴特沃斯低通滤波器函数求系统函数分子分母 , position* ,50,10,40

49、0,300）;%成一个窗口滋各窗口分为两行，第一部分用于绘制滤波系统幅频特性曲线滋会制滤波系统幅频特性曲线 悄示题滋制横纵坐标的标签癖二部分用于绘制滤波系统幅频特性曲线（db）燈制滤波系统幅频特性曲线（db）漪题淤制横纵坐标的标签遞波系统幅頻特性曲线000.511.522533.5数字频率(hrr数字頻率orr分析：这是设计的巴特沃斯低通滤波器的幅频特性曲线，观察该曲线可知该系统属于低通滤波器, 且符合设计指标的要求。2.5用数字滤波器进行滤波 调用filter函数溺用函数滤波sign aljvbo = filter(b z,az,sig nal_noise);figure( name*，滤波前后时域信号比较，卩osition* subplot(2,1,2);plot(tt,signal_lvbo);title('滤波后音频波形1 );xlabel( * 时间(s)' );ylabel(* 幅值')；subplot(2,1,1);plot( tt,signal_noise);title('滤