DSP课程设计-基于MATLAB的音乐信号处理和分析.doc

数字信号处理课程设计报告题目: 基于MATLAB的音乐信号处理和分析专业：测控技术与仪器年级：2011级姓名：学号：基于MATLAB的音乐信号处理和分析一、读取音乐信号并观察波形与频谱w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);sound(w,fs)figure;subplot(2,1,1);plot(wav);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);f=wwav/2*fs;subplot(2,1,2);plot(wwav,abs(fwav(1:lwav);1， 以二倍的抽样率听声音信号时，音乐播放的特别快，像被压缩了，播放的时间比原信号短。2， 以二分之一的抽样率听声音信号时，音乐播放的特别慢，像被拉长了，播放的时间比原信号长。3， 原信号频谱截止频率为0.5*pi 。二、（1）音乐信号的5倍减抽样w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);j=0;d1=5;for i=1:d1:lwav j=j+1; dwav1(j)=wav(i);endfigure;subplot(2,1,1);plot(dwav1);fdwav1=fft(dwav1);sound(dwav1,fs/d1)lwav1=round(length(fdwav1)/2);nwav=0:lwav1-1;wwav=nwav/(lwav1);f=wwav/2*fs;subplot(2,1,2);plot(wwav,abs(fdwav1(1:lwav1);%减抽样音乐频域（2）音乐信号的10倍减抽样（有混频）w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);j=0;d1=10;for i=1:d1:lwav j=j+1; dwav1(j)=wav(i);endfigure;subplot(2,1,1);plot(dwav1);fdwav1=fft(dwav1);sound(dwav1,fs/d1)lwav1=round(length(fdwav1)/2);nwav=0:lwav1-1;wwav=nwav/(lwav1);f=wwav/2*fs;subplot(2,1,2);plot(wwav,abs(fdwav1(1:lwav1);1， 原信号频谱截止频率为0.5*pi。2， 减抽样后的音乐信号听起来变得尖锐，有失真。3， 抽样率随着抽样间隔的增大而逐渐变小，声音越来越失真，音调变得急促，而尖锐，信号产生混叠3、 音乐信号的AM调制 （1）调制频率为pi/2w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);lam=length(wav);amw=pi/2;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;sound(amwav,fs)figure;subplot(2,1,1);plot(amwav);title(爸爸去哪儿调制后)famwav=fft(amwav);lwav=round(length(famwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);subplot(2,1,2);plot(wwav,abs(famwav(1:lwav);(2)调制频率为2*pi/3w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);lam=length(wav);amw=2*pi/3;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;sound(amwav,fs)figure;subplot(2,1,1);plot(amwav);title(爸爸去哪儿调制后)famwav=fft(amwav);lwav=round(length(famwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);subplot(2,1,2);plot(wwav,abs(famwav(1:lwav);1， 原信号的调制相当于频谱搬移，左移一个右移一个，当调制频率（余弦频率）小于0.4pi或大于0.6pi时就会产生混叠或丢失一部分信息。2， 当余弦点数取得少时，余弦频谱会产生泄漏。3， 当调制频率较高时（发生混叠），声音响度低，几乎只能听见兹兹的声音，信号几乎完全失真，当调制频率较低时（未发生混叠），声音很尖锐，响度较大，稍微能听出一点调子，但也有兹兹的声音。四、AM调制音乐信号的同步解调1.信号的解调w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);lam=length(wav);amw=pi/2;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;pmwav=amwav.*cwav;%音乐信号解调figure;subplot(2,1,1);plot(pmwav);title(解调时域)fpmwav=fft(pmwav);subplot(2,1,2);plot(wwav,abs(fpmwav(1:lwav);title(解调频域)sound(pmwav,fs);2.巴特沃斯滤波w,fs,bit=wavread(爸爸去哪儿.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);lam=length(wav);amw=pi/2;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;pmwav=amwav.*cwav;%音乐信号解调fpmwav=fft(pmwav);T=1;N,wc=buttord(0.05,0.5,1,15, s);B,A=butter(N,wc,s);%求模拟滤波器系统函数的系数C,D=bilinear(B,A,1/T);%求数字滤波器系统函数的系数W=0:0.001*pi:0.5*pi;H=freqs(B,A,W);%模拟滤波器的频率响应Hd=freqz(C,D,W);%数字滤波器频率响应figure;Hd_db=-20*log(abs(Hd(1)./(abs(Hd)+eps);subplot(3,1,1);plot(W/pi,abs(Hd);title(数字巴特沃斯滤波频响图)grid on;y=filter(C,D,pmwav);%对信号滤波fy=fft(y);subplot(3,1,2);plot(abs(y);title(滤波后音频)subplot(3,1,3);plot(wwav,abs(fy(1:lwav);title(滤波后频谱)sound(y,fs);3. 矩形窗滤波w,fs,bit=wavread(爸爸.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);lam=length(wav);amw=pi/2;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;pmwav=amwav.*cwav;%音乐信号解调fpmwav=fft(pmwav);N=39,wc=0.35*pi;%参数设置M=512;for n=0:N-1 if n=(N-1)/2 hd(n+1)=wc/pi; else hd(n+1)=sin(wc*(n-(N-1)/2)/(pi*(n-(N-1)/2); endend;%调用理想低通滤波器函数w2=boxcar(N);%设置矩形窗h2=hd.*w2;%求数字滤波器频率响应lx=round(length(h2)/2);wx1=(2/M)*0:M-1;fh2=fft(h2,M);db1=-20*log10(abs(fh2(1)./(abs(fh2)+eps);figure;subplot(4,1,1);plot(h2);title(矩形窗滤波特性)grid on;subplot(4,1,2);plot(wx1,abs(fh2);title(矩形窗频谱)grid on;y0=pmwav;yy2=conv(y0,h2);fyy2=fft(yy2);l2=round(length(fyy2)/2);wx2=(0:l2-1)/l2;subplot(4,1,3);plot(wx2,abs(yy2(1:l2);title(矩形窗滤波后音谱)grid on;subplot(4,1,4);plot(wx2,abs(fyy2(1:l2);title(矩形窗滤波后频谱)grid on;sound(yy2,fs);4、 布莱克曼窗滤波w,fs,bit=wavread(爸爸.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);lam=length(wav);amw=pi/2;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;pmwav=amwav.*cwav;%音乐信号解调fpmwav=fft(pmwav);N=43,wc=0.3*pi;%参数设置M=512;for n=0:N-1 if n=(N-1)/2 hd(n+1)=wc/pi; else hd(n+1)=sin(wc*(n-(N-1)/2)/(pi*(n-(N-1)/2); endend;%调用理想低通滤波器函数w2=blackman(N);%设置矩形窗h2=hd.*w2;%求数字滤波器频率响应fh2=fft(h2);lx=length(h2);wx1=(0:lx-1)*2/lx;db1=-20*log10(abs(fh2(1)./(abs(fh2)+eps);figure;subplot(4,1,1);plot(h2);title(布莱克曼窗滤波特性)grid on;subplot(4,1,2);plot(wx1,abs(fh2);title(布莱克曼窗频谱)grid on;y0=pmwav;yy2=conv(y0,h2);subplot(4,1,3);plot(yy2);title(布莱克曼窗滤波后音谱)grid on;fyy2=fft(yy2);l2=round(length(fyy2)/2);wx2=(0:l2-1)/l2;subplot(4,1,4);plot(wx2,abs(fyy2(1:l2);title(布莱克曼窗滤波后频谱)grid on;sound(yy2,fs);1， 解调后信号频谱在高频和低频处均有一部分，且成对称分布，需要滤掉高频才可大致还原原信号。2， 原信号的截止频率为0.5pi，使用数字巴特沃斯滤波器滤波器滤波参数通带截止频率0.5pi，阻带开始频率0.5pi，阻带衰减15db。滤波效果很好，基本还原了原信号。3， 使用窗函数滤波要根据过渡带宽算阶数N，选截止频率为0.5pi。4， 使用矩形窗滤波，矩形窗过渡带窄，但是阻带有波纹，高频部分有小部分未滤掉。5， 使用布莱克曼窗滤波，布莱克曼窗过渡带宽，但是阻带较好。6， 使用矩形窗和布莱克曼窗滤波，效果都行，基本都能还原原信号5、 音乐信号的滤波去噪1. 原始信号加入三余弦混合噪声w,fs,bit=wavread(爸爸.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);sf1=3000;sf2=5000;sf3=8000;l=length(wav);T=2/fs;t=0:T:(l-1)*T;s1=0.05*(cos(2*pi*sf2*t)+cos(2*pi*sf2*t)+cos(2*pi*sf3*t);%余弦噪声信号wav_s1=wav+s1;%加噪信号sound(wav_s1,fs)fwav_s1=fft(wav_s1);f_s1=fft(s1);figure;subplot(2,2,1);plot(s1);title(余弦噪声)subplot(2,2,2);plot(wav_s1);title(加噪信号)subplot(2,2,3);plot(wwav,abs(f_s1(1:lwav);subplot(2,2,4);plot(wwav,abs(fwav_s1(1:lwav);2. 用窗函数对混合噪声滤波w,fs,bit=wavread(爸爸.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);sf1=3000;sf2=5000;sf3=8000;l=length(wav);T=2/fs;t=0:T:(l-1)*T;s1=0.05*(cos(2*pi*sf2*t)+cos(2*pi*sf2*t)+cos(2*pi*sf3*t);%余弦噪声信号wav_s1=wav+s1;%加噪信号fwav_s1=fft(wav_s1);f_s1=fft(s1);N=45;wc=0.15*pi;%参数设置M=512;for n=0:N-1 if n=(N-1)/2 hd(n+1)=wc/pi; else hd(n+1)=sin(wc*(n-(N-1)/2)/(pi*(n-(N-1)/2); endend;%调用理想低通滤波器函数w1=hamming(N);%海明窗h1=hd.*w1;fh1=fft(h1);lx=length(h1);wx1=0:lx-1*2/lx;db1=-20*log10(abs(fh1(1)./(abs(fh1)+eps);figure;subplot(4,1,1);plot(h1);title(海明窗波形)grid on;subplot(4,1,2);plot(wx1,fh1);title(海明窗频谱)grid on;y0=wav_s1;yy1=conv(y0,h1);subplot(4,1,3);plot(yy1);title(海明窗滤波后音频)grid on;fyy1=fft(yy1);l2=round(length(fyy1)/2);wx2=(0:l2-1)/l2;subplot(4,1,4);plot(wx2,abs(fyy1(1:l2);title(海明窗滤波后频谱)grid on;sound(yy1,fs)3. 原始信号叠加随机白噪声并对其滤波w,fs,bit=wavread(爸爸.wav);wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);l=length(wav);s2=(rand(l,1)-0.5)/5;%用rand语句设置幅度为0.5的随机噪声wav_s2=wav+s2;fwav_s2=fft(wav_s2);f_s2=fft(s2);figure;subplot(2,2,1);plot(s2);title(随机白噪声)subplot(2,2,2);plot(wav_s2);title(加噪噪声)subplot(2,2,3);plot(wwav,abs(f_s2(1:lwav);subplot(2,2,4);plot(wwav,abs(fwav_s2(1:lwav);N=133;wc=0.15*pi;%参数设置M=512;fo