基于matlab 的音乐信号处理和分析

基于matlab的音乐信号处理和分析音乐信号的音谱和频谱的观察使用windows下的录音机录制一段音乐信号或者采用其他软件截取一段音乐信号（5s①使用d语句读取音乐信号获取抽样率（双列向量，需要分列处理；②输出音乐信号的波形和频谱，观察现象；③使用sound语句播放音乐信号，注意不同抽样率下的音调变化，解释现象程序：clearall;closeall;clc[y,fs,bit]=wavread('E:\music\11'); size(y) 看音乐信号是双列还是单列y1=y(:,1); 取单列fs=fs 获取音乐信号的抽样率N=length(y1); 音乐信号的长度Fy1=fft(y1,N); figure subplot(2,1,1);plot(y1);gridontitle('')xlabel('time/s')ylabel('Magnitude')gridontitle('音乐信号的频谱')xlabel('Frequency/pi')ylabel('Magnitude')sound(y1,fs) fs播放音乐信号sound(y1,2*fs) 2fs播放音乐信号播放音乐信号程序运行结果如下：音乐信号的抽样率fs=44100音乐信号的波形和频谱如下图：音乐信号的波形1uM

0.50-0.5-10 1 2 3 4time

5 6

8x1058000

音乐信号的频谱uM

60004000200000 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2抽样率较高时，听到的音乐音调较高，速度较快；抽样率较低时，听到的音乐音调低，速度慢。由此可知，音乐信号的音调和抽样率有关。音乐信号的抽取（减抽样）①（代表混叠与非混叠；②输出减抽样音乐信号的波形和频谱，观察现象，给出理论解释；③播放减抽样音乐信号，注意抽样率的改变，比较不同抽取间隔下的声音，解释现象。音乐信号的频率上限为3和0作为非混叠与混叠的抽样间隔进行抽样，程序如下：clearall;closeall;clc读取音乐信号y1=y(:,1); 取单列x1=y1(1:180000); d=3;j=0;fori=1:d:length(x1)j=j+1;x2(j)=x1(i);end 对音乐信号以3为间隔做减抽样N1=length(x1); 截断的音乐信号的长度N2=length(x2); 以3为间隔减抽样后音乐信号的长Fx1=fft(x1,N1); 对截断信号做快速傅里叶变换Fx2=fft(x2,N2); %对减抽样后的信号做快速傅里叶变figure subplot(2,1,1);plot(x1);gridontitle('')xlabel('time/s')ylabel('Magnitude')gridontitle('截断信号的频谱')xlabel('Frequency/pi')ylabel('Magnitude')figure 3为间隔减抽样后的信号的波形和频谱subplot(2,1,1);plot(x2);gridontitle('对截断信号以3为间隔减抽样后的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w2,abs(Fx2));gridontitle('对截断信号以3为间隔减抽样后的频谱')xlabel('Frequency/pi')ylabel('Magnitude')d=60;j=0;fori=1:d:length(x1)j=j+1;x3(j)=x1(i);end %60为间隔对截断信号进行减抽样N3=length(x3);Fx3=fft(x3,N3);w3=2/N3*[0:N3-1];figure %60为间隔减抽样后的音乐信号的波形和频谱subplot(2,1,1);plot(x3);gridontitle('对截断信号以60为间隔减抽样后的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w3,abs(Fx3));gridontitle('对截断信号以60为间隔减抽样后的频谱')xlabel('Frequency/pi')ylabel('Magnitude')sound(x1,fs) %播放截断信号%播放以3为间隔做减抽样后的音乐信号%播放以60为间隔做减抽样后的音乐信实验结果如图：截断后的波形1uM

0.50-0.5-10 2 4 6

time

10 12 14 16 18x104uM

300020001000

截断信号的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2对截断信号以3为间隔减抽样后的波形1uM

0.50-0.5-10 1 2 3time/s

4 5 6x104uM

1000500

对截断信号以3为间隔减抽样后的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2对截断信号以60为间隔减抽样后的波形1uM

0.50-0.5-10 500

1000

1500 2000 2500time/s

3000uM

对截断信号以60为间隔减抽样后的频谱60402000 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2使用两种抽样间隔进行实验：①以3为间隔进行减抽样，得到的音乐信号时域长度变为原来的1/3，频谱宽度变为原来的3倍，这是时域截断造成的频谱泄漏。由于原音乐信号的频率上限为0.2，故变为原来的3倍后仍不混叠，依然能够听到音乐；②以60为间隔进行减抽样，得到的音乐信号时域点数变为原来的1/60，相应的，频谱宽度变为原来的60倍，此时频谱已经混叠混叠，播放音乐信号时什么都听不到。AM调制①观察音乐信号的频率上限，选择适当调制频率对信号进行调制（率；②输出调制信号的波形和频谱，观察现象，给出理论解释；③播放调制音乐信号，注意不同调制频率下的声音，解释现象。0.2PI0.5pi0.95pi作为调制频率对音乐信号进行调制，程序如下：closeall;clearall;clc[y,fs,bit]=wavread('E:\music\22');y1=y(:,1);x1=y1(1:180000);余弦序列x4=x1.*c1'; %对信号进行调N4=length(x4); %调制后信号的长度Fx4=fft(x4,N4); %figure %subplot(2,1,1);plot(x4);gridontitle('对信号以频率pi/2调制后的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w4,abs(Fx4));gridontitle('')xlabel('Frequency/pi')ylabel('Magnitude')n1=1:180000;c1=cos(n1*pi*0.95);x5=x1.*c1';N5=length(x5);Fx5=fft(x5,N5);w5=2/N5*[0:N5-1];figuresubplot(2,1,1);plot(x5);gridontitle('对截断信号以频率pi*0.95调制后的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w5,abs(Fx5));gridontitle('对截断信号以频率pi*0.95调制后的频谱')xlabel('Frequency/pi')ylabel('Magnitude')sound(x1,fs)sound(x4,fs)sound(x5,fs)对信号以频率pi/2调制后的波形1uM

0.50-0.5-10 2 4 6

time/s

10 12 14 16 18x104uM

对信号以频率pi/2调制后的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2对截断信号以频率pi*0.95调制后的波形1uM

0.50-0.5-10 2 4 6

time/s

10 12 14 16 18x104uM

对截断信号以频率pi*0.95调制后的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2对音乐信号用两种频率进行调制：0.5pi处，不混叠，可以听到声音，但由于在高频处，音调很高，是刺耳的一段声音；调制频率为0.95pi时，调制后频谱搬移到0.95pi处，频谱混叠，不能听到声音。AM调制音乐信号的同步解调①计巴特沃斯滤波器完成同步解调；观察滤波器频率响应曲线；②FIR（形窗和布莱克曼窗，进行比较；③输出解调信号的波形和频谱，观察现象，给出理论解释；④播放解调音乐信号，比较不同滤波器下的声音，解释现象。设计巴特沃斯滤波器对解调信号进行滤波：clearall;closeall;clc[y,fs,bit]=wavread('E:\music\22');y1=y(:,1);x1=y1(1:180000);生成余弦序列x3=x1.*c1'; %调制x4=x3.*c1'; %解调N4=length(x4);w4=2/N4*[0:N4-1];figure %subplot(2,1,1);plot(x4);gridontitle('同步解调后音乐信号的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w4,abs(Fx4));gridontitle('同步解调后音乐信号的频谱')xlabel('Frequency/pi')ylabel('Magnitude')Wp=0.2;Ws=0.3;Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[B,A]=butter(N,Wc);[Hd,w]=freqz(B,A); %figureplot(w/pi,abs(Hd))gridontitle('数字巴特沃斯滤波器')xlabel('DigitalFrequency/pi')ylabel('Magnitude')y2=filter(B,A,x4); %N=length(y2);Fy2=fft(y1,N);w=2/N*[0:N-1];figuresubplot(2,1,1);plot(y2);gridontitle('巴特沃斯滤波后音乐信号的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w,abs(Fy2));gridontitle('巴特沃斯滤波后音乐信号的频谱')xlabel('Frequency/pi')ylabel('Magnitude')sound(x4,fs)sound(y2,fs)同步解调后音乐信号的波形1uM

0.50-0.5-10 2 4 6

time/s

10 12 14 16 18x104uM

同步解调后音乐信号的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2uM

1.41.210.80.60.40.2

数字巴特沃斯滤波器00 0.1 0.2 0.3 0.4 0.5 0.6DigitalFrequency/pi

0.7 0.8 0.9 1巴特沃斯滤波后音乐信号的波形0.4uM

0.20-0.2-0.40 2 4 6

time/s

10 12 14 16 18x104uM

300020001000

巴特沃斯滤波后音乐信号的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2用窗函数法设计矩形窗和布莱克曼窗对解调信号进行滤波：clearall;closeall;clc[y,fs,bit]=wavread('E:\music\22');size(y)y1=y(:,1);fs=fsx1=y1(1:90000);n1=1:90000;c1=cos(n1*pi/2);x3=x1.*c1';x4=x3.*c1';N4=length(x4);Fx4=fft(x4,N4);w4=2/N4*[0:N4-1];figuresubplot(2,1,1);plot(x4);gridontitle('同步解调后音乐信号的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w4,abs(Fx4));gridontitle('同步解调后音乐信号的频谱')xlabel('Frequency/pi')ylabel('Magnitude')N=33;wc=pi/4;hd=ideal(N,wc);w1=boxcar(N);w2=blackman(N);h1=hd.*w1';h2=hd.*w2';M=512;Fh1=fft(h1,M);Fh2=fft(h2,M);w=2/M*[0:M-1];figuresubplot(2,1,1);plot(w,abs(Fh1));gridontitle('矩形窗的频谱')xlabel('频率')ylabel('幅度')subplot(2,1,2);plot(w,abs(Fh2));gridon')xlabel('')Fy3=fft(y3,N4);figuregridontitle('')xlabel('time/s')ylabel('Magnitude')gridontitle('矩形窗滤波后的频谱')xlabel('Frequency/pi')ylabel('Magnitude')Fy4=fft(y4,N4);figuregridontitle('布莱克曼窗滤波后的波形')xlabel('time/s')ylabel('Magnitude')subplot(2,1,2);plot(w4,abs(Fy4));gridontitle('布莱克曼窗滤波后的频谱')xlabel('Frequency/pi')ylabel('Magnitude')1.51度幅0.5

矩形窗的频谱00 0.2 0.4

0.8

1 频率

1.4 1.6 1.8 21.51度幅0.5

布莱克曼窗的频谱00 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 2矩形窗滤波后的波形0.4uM

0.20-0.2-0.40 1 2 3 4 5time/s

6 7 8 9 10x104uM

1000500

矩形窗滤波后的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2布莱克曼窗滤波后的波形0.4uM

0.20-0.2-0.40 1 2 3 4 5time/s

6 7 8 9 10x104uM

1000500

布莱克曼窗滤波后的频谱00 0.2 0.4 0.6 0.8 1 Frequency/pi

1.4 1.6 1.8 2布莱克曼窗滤波效果更好，更平滑。但是听上去不同滤波器的效果差别不大。音乐信号的滤波去噪①原始信号叠加幅度为0.05，频率为3kHz,5kHz,8kHz的三余弦混合噪声，观察噪声频谱以及加噪后音乐信号的音频和频谱，并播放音乐，感受噪声对音乐信号的影响；②给原始音乐信号叠加幅度为5的随机白噪声（可用d语句产生，观察噪声频③根据步骤①、②观察到的频谱，选择合适指标设计滤波器进行滤波去噪，关注去噪后信号音谱和频谱，并播放音乐，解释现象。clearall;closeall;clc[y,fs,bit]=wavread('E:\music\22');y0=y(:,1);y1=y0(1:180000)%截取音乐信号N=length(y1);%取截断信号的长度n=0:N-1;x1=0.05*cos(n*2*pi*3000/fs)+0.05*cos(n*2*pi*5000/fs)+0.05*cos(n*2*pi*8000/fs);%三余弦噪声音乐信号和三余弦噪声叠加Fx1=fft(x1,N);Fy2=fft(y2,N);w1=2/N*[0:N-1];figuresubplot(3,1,1);plot(w1,abs(Fx1));gridontitle('三余弦混合噪声的频谱')xlabel('频率')ylabel('幅度')subplot(3,1,2);plot(w1,abs(Fy2));gridontitle('加三余弦混合噪声后音乐信号的频谱')xlabel('频率')ylabel('幅度')subplot(3,1,3);plot(y2);gridontitle('xlabel('')产生白噪声y3=y1+x2;%音乐信号和白噪声叠加Fx2=fft(x2,N);figuregridontitle('')xlabel('')subplot(3,1,2);plot(w1,abs(Fy3));gridontitle('加白噪声后音乐信号的频谱')xlabel('频率')ylabel('幅度')subplot(3,1,3);plot(y3);gridontitle('')xlabel('')Wp=0.09;Ws=0.11;Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[B,A]=butter(N,Wc);[Hd,w]=freqz(B,A);%设计巴特沃斯滤波器y4=filter(B,A,y2);%对加三余弦噪声的音乐信号进行滤波N2=length(y4);Fy4=fft(y4,N2);w2=2/N2*[0:N2-1];figuresubplot(2,1,1);plot(y4);gridontitle('对加三余弦混合噪声的音乐信号进行滤波之后的音谱')xlabel('时间')ylabel('幅度')subplot(2,1,2);plot(w2,abs(Fy4));gridontitle('xlabel('')Wp=0.15;Ws=0.17;Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[B,A]=butter(N,Wc);[Hd,w]=freqz(B,A);%设计巴特沃斯滤波器y5=filter(B,A,y3);%对加白噪声的音乐信号进行滤波N3=length(y5);Fy5=fft(y5,N3);w3=2/N3*[0:N3-1];figuresubplot(2,1,1);plot(y5);gridontitle('对加白噪声的音乐信号进行滤波之后的音谱')xlabel('时间')ylabel('幅度')subplot(2,1,2);plot(w2,abs(Fy5));gridon')xlabel('')sound(y2,fs)sound(y3,fs)sound(y4,fs)sound(y5,fs)5000度幅0

三余弦混合噪声的频谱0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2频率5000度幅0

加三余弦混合噪声后音乐信号的频谱0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2频率加三余弦混合噪声后音乐信号的音谱1度 0幅-10 2 4 6 8时间

10 12 14 16 18x104400度幅

白噪声的频谱00 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 24000度幅

加白噪声后音乐信号的频谱00 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 2加白噪声后音乐信号的音谱1度 0幅-10 2 4 6 8时间

10 12 14 16 18x1041度 幅

对加三余弦混合噪声的音乐信号进行滤波之后的音谱-0.5-10 2 4 6 8时间

10 12 14 16 18x10430002000度幅1000

对加三余弦混合噪声的音乐信号进行滤波之后的频谱00 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 2对加白噪声的音乐信号进行滤波之后的音谱10.5度 0幅-0.5-10 2 4 6 8时间

10 12 14 16 18x10430002000度幅1000

对加白噪声的音乐信号进行滤波之后的频谱00 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 2对于加三余弦混合噪声的音乐信号，播放的过程中，背后一直有一个音调很高的“滴”的声音，这是由于余弦信号频率较高；对于加白噪声的音乐信号，播放过程中，背后一直有“丝啦丝啦”的干扰噪声，就像信号不好时听收音机的声音；根据滤波后的音乐播放效果可知，三余弦混合噪声可以完全滤去，这是由于余弦信号不在音乐信号所在的频率处；白噪声却不能完全滤去，这是由于白噪声在低频处也有分布，音乐信号也存在低频处，因此低频处的白噪声无法滤除，所以听到的音乐信号仍然会有干扰的声音，但是效果比滤波前要好。音乐信号的幅频滤波及相频分析①（可自行选择不同的截止频率影响；②（可自行选择不同的截止频率影响；③选取两端不同的音乐信号，分别将其幅度谱与相位谱交叉组合构成信的音乐信号，播放比较组合后的音乐与原始音乐，感受相频信息对音乐信号的影响。clearall;closeall;clc[y,fs,bit]=wavread('E:\music\22');y1=y(:,1);Wp=0.06;Ws=0.09;Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[B,A]=butter(N,Wc);[Hd,w]=freqz(B,A);%设计巴特沃斯低通滤波器y2=filter(B,A,y1);%滤除音乐信号的高频信息N1=length(y1);N2=length(y2);Fy1=fft(y1,N1);Fy2=fft(y2,N2);w1=2/N1*[0:N1-1];w2=2/N2*[0:N2-1];figure%绘制低通滤波前后音乐信号的幅頻特性subplot(2,1,1);plot(w1,abs(Fy1));gridontitle('低通滤波前音乐信号的频谱')xlabel('频率')ylabel('幅度')subplot(2,1,2);plot(w2,abs(Fy2));gridontitle('')xlabel('')Wp=0.06;Ws=0.09;Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[B,A]=butter(N,Wc,'high');[Hd,w]=freqz(B,A);%设计巴特沃斯高通滤波器y3=filter(B,A,y1);%滤除音乐信号的低频信息N3=length(y3);w3=2/N3*[0:N3-1];Fy3=fft(y3,N3);figure%绘制高通滤波前后音乐信号的幅頻特性subplot(2,1,1);plot(w1,abs(Fy1));gridontitle('高通滤波前音乐信号的频谱')xlabel('频率')ylabel('幅度')subplot(2,1,2);plot(w3,abs(Fy3));gridontitle('')xlabel('')sound(y1,fs)sound(y2,fs)sound(y3,fs)低通滤波前音乐信号的频谱80006000度4000幅200000 0.2

0.6

1 1.2 频率

1.6 1.8 280006000度幅

低通滤波后音乐信号的频谱200000 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 280006000度幅

高通滤波前音乐信号的频谱200000 0.2

0.6

1 1.2 频率

1.6 1.8 215001000度幅500

高通滤波后音乐信号的频谱00 0.2 0.4 0.6 0.8 1频率

1.2 1.4 1.6 1.8 2低通滤波后的音乐信号音调较原音乐信号低，高通滤波后的音乐信号音调较原音乐信号高。clearall;closeall;clc[x1,fs1,bit1]=wavread('E:\music\44');[x2,fs2,bit2]=wavread('E:\music\22');y11=x1(:,1);y22=x2(:,1);N1=length(y1);Fy1=fft(y1,N1);N2=length(y2);Fy2=fft(y2,N2);figure %画两信号的幅度谱gridontitle('44')xlabel('Frequency/Hz')ylabel('Magnitude')gridontitle('音乐信号22的幅度谱')xlabel('Frequency/Hz')ylabel('Magnitude')subplot(2,2,3);plot(w1,angle(Fy1));gridon %画两信号的相位谱title('音乐信号44的相位谱')xlabel('频率')ylabel('相位')subplot(2,2,4);plot(w2,angle(Fy2));gridontitle('22')xlabel('')Fy3=abs(Fy1).*exp(1i*a