matlab处理音乐数字信号.doc

数字信号处理学院：物理工程学院 姓名：金义飞 1、 读取音乐信号并观察波形与频谱程序w,fs,bit=wavread(你是我的眼.wav); %音频信号w。采样率fs 精度bit wav=(w(:,1); %取第一列信号 sound(w,fs) %听取抽样后的wav音乐sound(w,2*fs) %听取抽样后的wav音乐sound(w,fs/2) %听取抽样后的wav音乐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; %x轴坐标 实际频率 subplot(2,1,2);plot(wwav,abs(fwav(1:lwav); %画出音乐频域图图像现象说明：从音谱即时域谱可以看出音乐信号随时间的幅度变化；从频谱可看到音乐信号主要分布在低频段，高频成分较少，在0.5pi 以后几乎无音乐信号2、 音乐信号的抽取（减抽样）三倍减抽样1，程序clear all;close all;clc %清除存储空间，关闭所有窗口，清除命令w,fs,bit=wavread(你是我的眼.wav)wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);j=0;d1=3; for i=1:d1:lwav; j=j+1; dwav1(j)=wav(i);end %本段每隔三个点取一个点，相当于三倍减抽样sound(dwav1,fs/d1);figure;subplot(2,1,1);plot(dwav1); fwav=fft(dwav1);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);f=wwav/fs;subplot(2,1,2);plot(wwav,abs(fwav(1:lwav);2图像20倍减抽样1，程序clear all;close all;clcw,fs,bit=wavread(你是我的眼.wav)wav=(w(:,1);fwav=fft(wav);lwav=round(length(fwav)/2);j=0;d1=20; for i=1:d1:lwav; j=j+1; dwav1(j)=wav(i);endsound(dwav1,fs/d1);figure;subplot(2,1,1);plot(dwav1); fwav=fft(dwav1);lwav=round(length(fwav)/2);nwav=0:lwav-1;wwav=nwav/(lwav);f=wwav/fs;subplot(2,1,2);plot(wwav,abs(fwav(1:lwav); 分别选取抽样间隔为3倍频和20倍频进行减抽样，后者产生混叠，前者不混叠，当选取3倍频减抽样时，频率样间隔变大但未超过最小抽样间隔，尽管频谱有所延伸，但并不产生混叠，即减抽样后频谱上限未超过折叠频率；而选取20倍频减抽样时，由于时域抽样间隔变大超过最小抽样间隔，频谱延伸且产生混叠，即频率上限超过折叠频率，与3倍减抽样相比声音信号中多了很多杂音，这是由于抽样间隔的增大使信号频谱向高频搬移，且间隔越大，搬移越厉害。 3、 音乐信号的AM调制1、 调制频率wc=0.7*pi程序w,fs,bit=wavread(你是我的眼.wav);wav=(w(:,1);lam=length(wav);fwav=fft(wav);lwav=length(fwav);nwav=0:lwav-1;amw=0.7*pi;cwav=cos(amw*0:lam-1);amwav=wav.*cwav;sound(amwav,fs);figure;subplot(2,1,1);plot(amwav);title(你是我的眼时域调制)famwav=fft(amwav);subplot(2,1,2);wwav=nwav/lwav;plot(2*wwav,abs(famwav(1:lwav);title(你是我的眼频域调制)图像2、 调制频率wc=pi/2程序w,fs,bit=wavread(你是我的眼.wav);wav=(w(:,1);lam=length(wav);fwav=fft(wav);lwav=length(fwav);nwav=0:lwav-1;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);subplot(2,1,2);wwav=nwav/lwav;plot(2*wwav,abs(famwav(1:lwav);title(你是我的眼频域调制)图像声音原信号频率上限大约为0.4pi,选取高频调制频率为0.7pi,低频调制频率为Pi/2; 高频调制后，信号频谱从原点处搬移至0.7pi处，且频谱左右各搬移0.7pi,由于频率上限为0.4pi,故搬移后信号最高频率超过折叠频率，从而产生了混叠，将会使信号失真;低频调制后，信号频谱从原点搬移至pi/20处，且频谱左右各搬移Pi/2,但搬移后信号频率上限并未超出折叠频率，故信号未发生混叠；播放调制后的声音信号，可以听出高频调制后声音信号有刺耳的噪音，而低频调制仍能清楚的听出其音调。因为高频调制产生混叠而低频没有产生混叠。 四、AM调制音乐信号的同步解调1、巴特沃斯滤波程序clear all;close all;clcw,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);figure;subplot(2,1,1);plot(pmwav);title(解调时域)fpmwav=fft(pmwav);subplot(2,1,2);plot(wwav,abs(fpmwav(1:lwav);title(解调频域)T=1;wp=2/T*tan(0.3*pi/2);ws=2/T*tan(0.4*pi/2);Rp=1;Rs=15; N,wc=buttord(wp,ws,Rp,Rs,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);%数字滤波器频率响应figuresubplot(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);图像2、矩形窗滤波clear all;close all;clcw,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;hd=ideal(N,wc);w2=boxcar(N);%设置矩形窗h2=hd.*w2;%求滤波器频率响应fh2=fft(h2,M);lx=length(fh2);wx1=(2/M)*0:M-1;db1=-20*log10(abs(fh2(1)./(abs(fh2)+eps);figure;subplot(3,1,1);plot(wx1,db1);title(矩形窗滤波特性);grid on;y0=pmwav;yy2=conv(y0,h2);subplot(3,1,2);plot(yy2);title(矩形窗滤波后音谱);grid on;fyy2=fft(yy2);l2=round(length(fyy2)/2);wx2=(0:l2-1)/l2;subplot(3,1,3);plot(wx2,abs(fyy2(1:l2);title(矩形窗滤波后频谱);grid on;sound(yy2,fs);定义函数文件function hd=ideal(N,wc)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图像4、布莱克曼窗滤波clear all;close all;clcw,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=75;wc=0.15*pi;hd=ideal(N,wc);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(3,1,1);plot(wx1,db1);title(blackman窗滤波特性);grid on;y0=pmwav;yy2=conv(y0,h2);subplot(3,1,2);plot(yy2);title(blackman窗滤波后音谱);grid on;fyy2=fft(yy2);l2=round(length(fyy2)/2);wx2=(0:l2-1)/l2;subplot(3,1,3);plot(wx2,abs(fyy2(1:l2);title(blackman窗滤波后频谱);grid on;sound(yy2,fs);函数文件function hd=ideal(N,wc)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图像分析：巴特沃斯滤波器的特点是通带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。五、音乐信号的滤波去噪1、原始信号叠加余弦混合噪声clear all;close all;clcw,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*sf1*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,3);plot(wwav,abs(f_s1(1:lwav);title(余弦噪声频域);subplot(2,2,2);plot(wav_s1);title(加噪信号时域);subplot(2,2,4);plot(wwav,abs(fwav_s1(1:lwav);title(加噪信号频域);用巴特沃斯滤波clear all;close all;clcw,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*sf1*t)+cos(2*pi*sf2*t)+cos(2*pi*sf3*t);wav_s1=wav+s1;T=1;wp=2/T*tan(0.15*pi/2);ws=2/T*tan(0.2*pi/2);Rp=1;Rs=10; N,wc=buttord(wp,ws,Rp,Rs,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);%数字滤波器频率响应figuresubplot(3,1,1);plot(W/pi,abs(Hd);title(数字巴特沃斯滤波特图)grid on;y=filter(C,D,wav_s1);%对信号滤波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、原始信号叠加随机白噪声并对其滤波。clear all;close all;clcw,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(1,l)-0.5)/5;%用rand语句设置幅度为0.5的随机白噪声f_s2=fft(s2);wav_s2=wav+s2;fwav_s2=fft(wav_s2);figure;subplot(2,2,1);plot(s2);title(随机白噪声时域);subplot(2,2,3);plot(wwav,abs(f_s2(1:lwav);title(随机白噪声频域);subplot(2,2,2);plot(wav_s2);title(加噪噪声时域);subplot(2,2,4);plot(wwav,abs(fwav_s2(1:lwav);title(加噪噪声频域);图片程序clear all;close all;clcw,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(1,l)-0.5)/5;%用rand语句设置幅度为0.5的随机白噪声f_s2=fft(s2);wav_s2=wav+s2;fwav_s2=fft(wav_s2);N=39;wc=0.35*pi;%参数设置M=512;hd=ideal(N,wc);w2=boxcar(N);%设置矩形窗h2=hd.*w2;%求滤波器频率响应fh2=fft(h2,M);lx=length(fh2);wx1=(2/M)*0:M-1;db1=-20*log10(abs(fh2(1)./(abs(fh2)+eps);figure;subplot(3,1,1);plot(wx1,db1);title(矩形窗滤波特性);grid