实验三音频信号分析与处理-3文档

长春理工大学国家级电工电子实验教学示范中心学生实验报告——学年第学期实验课程实验地点学院专业学号姓名

实验项目实验时间实验台号预习成绩报告成绩一、实验目的掌握音频信号的采集方法掌握运用Matlab实现对音频信号的频谱分析掌握运用Matlab设计RC滤波系统掌握运用Matlab实现对加噪声后的音频信号的进行滤波处理的方法二、实验性质设计性实验三、实验任务音频信号的采集音频信号的采集可以通过Windows自带的录音机也可以用专用的录制软件录制一段音频信号（尽量保证无噪音、干扰小），也可以直接复制一段音频信号，但必须保证音频信号保存为.wav的文件。音频信号的时域、频域分析运用Matlab软件实现对音频信号的打开操作及信号的时域分析和信号的频域分析，并画出相应的图形（要求图形有标题），并打印在实验报告中（注意：把打印好的图形剪裁下来，粘贴到实验报告纸上）。引入干扰信号在原有的音频信号基础上，叠加一个频率为100KHz的正弦波干扰信号（幅度自定，可根据音频信号的情况而定）。滤波系统的设计运用Matlab实现RC滤波系统，要求加入干扰信号的音频信号经过RC滤波系统后，能够滤除100KHz的干扰信号，同时保留原有的音频信号，绘制出RC滤波系统的冲激响应波形并分析其频谱。四、实验设计分析原信号分析：首先利用MATLAB中函数audioread即[y,Fs]=audioread('挪威的森林.wav');读取音频文件，其中Fs为采样频率且默认为44100hz，用SampleRate=get(hAudio,'SampleRate');%获取音频文件的采样率KHz，利用myfft可实现傅里叶变变换求出原信号的频谱，从频谱图可以看出原信号的频谱分布主要在(0,7.5*10^3)Hz并用plot函数输出用play函数播放原音频。引入噪声：根据叠加原理可引入100khz的噪声并用sin（2π*10^5）表示，同样可以myfft可实现傅里叶变变换求出加入噪声后信号的频谱,可知加入噪声后其频谱中引入了频率约为12kHz的频率分量，是需要滤除的部分。RC滤波：根据以上分析需滤掉的频率分量高于原信号的频率分量故此处采用二阶有源低通滤波器，二阶有源低通滤波器公式H=A0*Wc^2./((j*w).^2+Wc/Q*(j*w)+Wc^2);由于原信号的频谱分布主要在(0,7.5*10^3)Hz，而引入噪声的频谱频率约为12kHz根据二阶有源低通滤波器的原理及幅频响应曲线当等效品质因数Q较小滤波效果更好故Q选为0.5，经试验当选取特征角频率为4000时其滤波效果较好且原信号频谱失真较小，则根据Q=1/(3-A0)得放大倍数A0=1;开始（二）流程图开始添加文件读取音频信号添加文件读取音频信号原信号时域及频域分析时域图原信号时域及频域分析时域图播放原信号播放原信号频域图频域图100kHz噪声100kHz噪声加入噪声后信号的频谱分析加入噪声后信号的频谱分析播放加入噪声后的文件频域图播放加入噪声后的文件频域图滤波器的幅频响应曲线滤波器的幅频响应曲线通过RC滤波器后的频谱分析通过RC滤波器后的频谱分析滤波后的频谱图滤波后的频谱图播放滤波后的音频滤除噪声后傅里叶变换得原信号播放滤波后的音频滤除噪声后傅里叶变换得原信号时域图时域图结束结束（三）结果及分析从原音频的时域及频域图可以看出原信号的频谱分布主要在(0,7.5*10^3)Hz引入噪声后其频谱中引入了频率约为12kHz的频率分量，是需要滤除的部分。由设计的二阶有源低通滤波器的幅频响应曲线可知其对5khz以上的频率有较好的滤除作用。由滤波以后的频谱可以看出其较好的滤除了噪声而保留了原信号。滤波后的频谱的傅里叶变换得到的时域波形与原信号的时域波形几乎一致，说明滤波效果较好。实验程序%音频信号分析与处理%%打开和读取音频文件clearall;%清除工作区缓存[y,Fs]=audioread('sh.wav');%读取音频文件VoiceWav=y(300000:400000,1);%截取音频中的一段波形cleary;%清除缓存hAudio=audioplayer(VoiceWav,Fs);%将音频文件载入audioplayerSampleRate=get(hAudio,'SampleRate');%获取音频文件的采样率KHzT=1/SampleRate;%计算每个点的时间，即采样周期SampLen=size(VoiceWav,1);%单声道采样长度%%绘制时域分析图hFig1=figure('Units','normalized','Position',[00.050.490.85]);t=T:T:(SampLen*T);subplot(2,1,1);%绘制音频波形plot(t,VoiceWav);%绘制波形title('音频时域波形图');axis([0,2.3,-0.5,0.5]);xlabel('时间(s)');ylabel('幅值(V)');%显示标题%%傅里叶变换subplot(2,1,2);%绘制波形myfft(VoiceWav,SampleRate,'plot');%傅里叶变换title('单声道频谱振幅');%显示标题xlabel('Frequency(Hz)');1ylabel('|Y(f)|');play(hAudio);%播放添加噪声前的声音pause(3);%%引入100KHz的噪声干扰t=(0:SampLen-1)*T;noise=0.06*sin(2*pi*100000*t);%噪声频率100Khz，幅值-0.06V到+0.06VhFig2=figure('Units','normalized','Position',[0.50.050.50.85]);subplot(2,1,1);%绘制波形plot(t(1:1000),noise(1:1000));title('100KHz噪声信号');%显示标题noiseVoice=VoiceWav+noise';%将噪声加到声音里面hAudio=audioplayer(noiseVoice,Fs);%将音频文件载入audioplayersubplot(2,1,2);%绘制波形[fftNoiseVoice,f]=myfft(noiseVoice,SampleRate,'plot');title('音乐和噪声频谱');%显示标题play(hAudio);%播放添加噪声后的声音pause(3);%%设计RC滤波系统（二阶有源低通滤波器）w=f;Wc=3000;%wc=1/(RC)，特征角频率A0=1;%A0=AVF<3Q=1/(3-A0);%品质因素H=A0*Wc^2./((j*w).^2+Wc/Q*(j*w)+Wc^2);%二阶有源低通滤波器公式hFig3=figure('Units','normalized','Position',[00.050.490.85]);subplot(2,1,1);%绘制波形plot(w,H);%显示标题title('二阶有源低通滤波器');xlabel('频率');ylabel('w/Wc');AfterFilter=fftNoiseVoice.*H';%滤波%AfterFilter=fftNoiseVoice;%%傅里叶逆变换subplot(2,1,2);%绘制波形NFFT=2^nextpow2(SampLen);%根据采样求傅里叶变换的点f=SampleRate/2*linspace(0,1,NFFT/2+1);%计算频率显示范围plot(f,2*abs(AfterFilter(1:NFFT/2+1)));%绘制频域分析图title('滤波以后的频谱');xlabel('频率');ylabel('w/Wc');clearYH;FilterVoice=fftshift(ifft(AfterFilter));%傅里叶逆变换FilterVoice=fftshift(FilterVoice);FilterVoice=((FilterVoice-max(max(FilterVoice)))/(max(max(FilterVoice))...-min(min(FilterVoice))))+0.5;%归一化clearAfterFilterw;%清除缓存hFig4=figure('Units','normalized','Position',[0.50.050.50.85]);subplot(2,1,1);%绘制右声道波形plot(t,FilterVoice(1:size(t,2),1));title('傅里叶逆变换图');xlabel('时间(ms)');ylabel('幅值(V)');%显示标题hAudio=audioplayer(FilterVoice,Fs);%将音频文件载入audioplayerplay(hAudio);%播放滤掉噪声后的声音%%结束傅里叶变换子函数function[outFFT,Freq]=myfft(varargin)%输入参数格式：%1.需要FFT变换的向量%2.采样率SampleRate%3.是否绘图，绘图‘plot‘，不绘图则不传递该参数%4.单边显示：'half'，全部显示：'full'%输出参数格式：%1.转换完成的向量%2.FFT频率范围%%输出参数判断switchnargincase0|1,error('Lessargumentin!');case2,FFTVector=varargin{1};SampleRate=varargin{2};isplot=0;case3,FFTVector=varargin{1};SampleRate=varargin{2};isplot=varargin{3};plotmode='half';case4,FFTVector=varargin{1};SampleRate=varargin{2};isplot=varargin{3};plotmode=varargin{4};otherwiseerror('Somanyargumentsin!');end%%FFT变换SampLen=size(FFTVector,1);NFFT=2^nextpow2(SampLen);Y=fft(FFTVector,NFFT)/SampLen;f=SampleRate*linspace(0,1,NFFT);%%判断输出参数if(nargout==1)outFFT=Y;elseif(nargout==2)outFFT=Y;Freq=f;end%获取采样点%根据采样求傅里叶变换的点%傅里叶变换%计算频率显示范围%%判断绘图if(strcmp(isp