数字信号处理课程设计基于 matlab 的音乐信号处理和分析

设计题目：基于 MATLAB 的音乐信号处理和分析

一、课程设计的目的

    本课程设计通过对音乐信号的采样、抽取、调制、解调等多种处理过程的理论分析和MATLAB实现，使学生进一步巩固数字信号处理的基本概念、理论、分析方法和实现方法；使学生掌握的基本理论和分析方法知识得到进一步扩展；使学生能有效地将理论和实际紧密结合；增强学生软件编程实现能力和解决实际问题的能力。

二、课程设计基本要求

1 学会 MATLAB 的使用， 掌握MATLAB的基本编程语句。

2 掌握在 Windows 环境下音乐信号采集的方法。

3 掌握数字信号处理的基本概念、基本理论和基本方法。

4 掌握 MATLAB 设计 FIR 和 IIR 数字滤波器的方法。 

5 掌握使用MATLAB处理数字信号、进行频谱分析、设计数字滤波器的编程方法。

三、课程设计内容

1、音乐信号的音谱和频谱观察

使用windows下的录音机录制一段音乐信号或采用其它软件截取一段音乐信号（要求：时间不超过5s、文件格式为wav文件）

① 使用wavread语句读取音乐信号，获取抽样率；（注意：读取的信号是双声道信号，即为双列向量，需要分列处理）；

② 输出音乐信号的波形和频谱，观察现象；

③ 使用sound语句播放音乐信号，注意不同抽样率下的音调变化，解释现象。

Wavread格式说明：

[w,fs,b]=wavread(‘语音信号’)，采样值放在向量w中，fs表示采样频率（hz），b表示采样位数。

上机程序：

[y,fs,bit]=wavread('I do片段')%读取音乐片段，fs是采样率

size(y)%求矩阵的行数和列数

y1=y( : ,1);%对信号进行分列处理

n1=length(y1);%取y的长度

t1=(0:n1-1)/fs;%设置波形图横坐标

figure

subplot(2,1,1);

plot(t1,y1);      %画出时域波形图

ylabel('幅值');

xlabel('时间（s）');

title('信号波形');

subplot(2,1,2);

Y1=fft(y1);

w1=2/n1*(0:n1-1);%设置角频率

plot(w1,abs(Y1));%画频谱图

title('信号频谱');

xlabel('数字角频率');

ylabel('幅度');

grid on;

sound(y,fs); 

实验结果：

1、通过观察频谱知，选取音乐信号的频谱集中在0~0.7*pi之间，抽样点数fs=44100;

2、当采样频率问原来0.5（0.5*fs）倍时：音乐片段音调变得非常低沉，无法辨认原声，播放时间变长；抽样频率减小，抽样点数不变时，其分辨力增大，记录长度变长，声音失真。

3、当采样频率问原来2(2*fs)倍时：音乐片段音调变得尖而细，语速变快，播放时间变短；抽样频率增加，抽样点数不变时，其分辨力下降，记录长度变短，声音失真。

2、音乐信号的抽取（减抽样）

① 观察音乐信号频率上限，选择适当的抽取间隔对信号进行减抽样（给出两种抽取间隔，代表混叠与非混叠）；

② 输出减抽样音乐信号的波形和频谱，观察现象，给出理论解释； 

③ 播放减抽样音乐信号，注意抽样率的改变，比较不同抽取间隔下的声音，解释现象。

上机程序：

[y,fs,bit]=wavread('I do片段')

y1=y( : ,1);

n1=length(y1);

tn1=(0:n1-1)/fs;

figure 

subplot(2,1,1);

plot(tn1,y1);    

ylabel('幅度');

xlabel('时间（s）');

title('原信号波形');

wn1=2/n1*[0:n1-1];

Y1=fft(y1); 

subplot(2,1,2);

plot(wn1,abs(Y1));

title('原信号频谱');

xlabel('数字角频率w');

ylabel('幅度');

grid on;

D=2;%设置抽样间隔 

y2=y1(1:D:n1);%减抽样

n2=length(y2);%减抽样后信号长度

t2=(0:n2-1)/fs;%设置横坐标

figure 

subplot(2,1,1);

plot(t2,y2);     %绘制减抽样信号波形图

ylabel('幅度');

xlabel('时间（s）');

title('2:1减抽样信号波形');

Y2=fft(y2); %对y2进行n2点fft谱分析

w2=2/n2*[0:n2-1];

subplot(2,1,2);

plot(w2,abs(Y2));%绘制减抽样信号频谱图

title('2:1减抽样信号频谱');

xlabel('数字角频率w');

ylabel('幅度');

grid on;

sound(y2,fs/D);

实验结果与分析：

1、程序中指标D表示抽样间隔，其值越大，相邻两抽样点之间的距离越远，抽样后漏掉的信息越多，相应的时域信号长度越短；

2、 抽样间隔D=1.1时的信号波形及频谱图，抽样频率大于信号最高频率的两倍，

满足抽样定理，不会发生混叠。抽样间隔D越大，抽样率fs越小，抽样后时域信号长度越短

3 、抽样间隔D=2时的信号波形及频谱图  

抽样间隔D=2的信号波形及频谱图，抽样频率小于信号最高频率的两倍，即<2，不满足抽样定理，其频谱图发生混叠，且D越大，混叠越严重，高频成分增加越多，音乐片段音调听起来很沙哑， 音调变得很高。

4、抽样间隔D=5的信号波形及频谱图，抽样频率小于信号最高频率的两倍，即<2，不满足抽样定理，其频谱图发生混叠，音乐片段音调听起来很尖锐。

3、音乐信号的AM调制

① 观察音乐信号频率上限，选择适当调制频率对信号进行调制（给出高、低两种调制频率）；

② 输出调制信号的波形和频谱，观察现象，给出理论解释； 

③ 播放调制音乐信号，注意不同调制频率下的声音，解释现象。

上机程序：

[y,fs,bit]=wavread('I do片段')

y1=y( : ,1);

n1=length(y1);

n3=0:(n1-1);

b1=cos(0.75*pi*n3);%设置调制信号  

c1=b1'.*y1;%对原信号进行调制

lc1=length(c1);

t=(0:lc1-1)/fs;

figure                   %用载波对信号进行调制，并对其做fft变换

subplot(2,1,1)          %获取频谱，从图中可以观察到，调制后的

plot(t,c1);                %信号频谱发生搬移

xlabel('时间(s)');

ylabel('幅度');

title('调制后信号');

w1=2/lc1*[0:lc1-1];%设置角频率W

C1=fft(c1);

subplot(2,1,2)

plot(w1,abs(C1));

xlabel('数字角频率w');

ylabel('幅度');

title('调制后信号的频谱(高频率调制)');

grid on;

%sound(c1,fs); 

实验结果与分析：

信号的调制过程就是将信号频谱搬移到任何所需的较高频率范围。调制的实质是把各种信号的频谱搬移，使它们互不重叠的占据不同的频率范围，也即信号分别托付于不同频率的载波上。具体的调制原理推导在此不再叙述，仅将结论列出：f(t)=g(t)*cos(w0*t)。由此将信号g(t)的频谱搬移到（2*n+1）w附近，同时音乐信号的频谱幅度变为原来的1/2。

如果信号的最高频谱wh超过了ws/2，则各周期延括分量产生频谱的交叠，称为是频谱的混叠现象。根据奈奎斯特定律可知，若希望频谱不会发生混叠，则fs>=fh。

 调制后信号波形及频谱（低频率调制）b1=cos(0.25*pi*n3);

 调制后信号波形及频谱（高频率调制）b1=cos(0.75*pi*n3);

1、由频谱图知信号频率上限约在0.7*pi处，取高调制频率为0.75*pi，取高调制频率为0.25*pi；观察调制后的频谱知，高频率调制后频谱混叠，而低频率调制后频谱未发生混叠；

2、用载波对信号进行调制，其实质过程为在时域用载波函数b1=cos(0.25*pi*n3)乘上原信号函数，在频域既表现为信号频谱的搬移，且载波函数角频率越大，调制后信号频谱搬移越大；

3、调制后信号高频成分增加，且调制频率越高，高频成分越大，声音变得越低沉，失真程度越大。

4、AM调制音乐信号的同步解调

① 设计巴特沃斯IIR滤波器完成同步解调；观察滤波器频率响应曲线；

② 用窗函数法设计FIR滤波器完成同步解调，观察滤波器频率响应曲线；（要求：分别使用矩形窗和布莱克曼窗，进行比较）； 

③ 输出解调信号的波形和频谱，观察现象，给出理论解释； 

④ 播放解调音乐信号，比较不同滤波器下的声音，解释现象。

上机程序：

clear all;cla;close all

[a,fs,bit]=wavread('I do片段');

y1=a( : ,1);%去单列数据进行分析

f1=fft(y1);

n=length(f1);

tn=(0:n-1)/fs;

w=2/n*[0:n-1];

%sound(y1,fs);

figure(1)

subplot(2,2,1);plot(tn,y1);%绘制音乐信号时域波形图

grid on;

title('原信号音频');

xlabel('时间');

ylabel('幅度');

subplot(2,2,2);plot(w,abs(f1));%绘制音乐信号频谱图

grid on;

title('原信号频谱');

xlabel('频率/pi');

ylabel('幅度');

t=[0:n-1];

y2=cos(pi*1/2*t);%载波函数

y3=y1.*y2';%信号调制

ty3=(0:length(y3)-1)/fs;

subplot(2,2,3);plot(ty3,y3);%绘制调制后信号波形图

grid on;

title('AM调制音频信号');

xlabel('时间');

ylabel('幅度');

f3=fft(y3);

n2=length(f3);

w2=2/n2*[0:n2-1];

subplot(2,2,4);plot(w2,abs(f3));%绘制调制后信号频谱图

grid on;

title('AM调制频谱');

xlabel('频率/pi');

ylabel('幅度');

%解调后信号

y4=y3.*y2'%信号解调

ly4=length(y4);

ty4=(0:(ly4-1))/fs;

wy4=2/ly4*[0:ly4-1];

Y4=fft(y4);

figure(2)

subplot(3,1,1);plot(ty4,y4);

grid on;

title('解调后的波形');

xlabel('t');

ylabel('幅度');

subplot(3,1,2);plot(wy4,abs(Y4));

title('解调后波形频谱');

xlabel('w/pi');

ylabel('幅度');

grid on;

%设计巴特沃斯滤波器进行滤波去噪

[N1,wc1]=buttord(0.05,0.50,1,15);%确定低通滤波器的阶数和截止频率；

[b,a]=butter(N1,wc1);           %确定低通滤波器分子分母系数

[H,W]=freqz(b,a);

subplot(3,1,3);

plot(W,abs(H));

title('IIR滤波器频谱');

xlabel('w/pi');

ylabel('幅度');

grid on;

m=filter(b,a,y4);

sound(m,fs);

lm=length(m);%滤波后信号长度

tm=(0:lm-1)/fs;%设置横坐标

wm=2/lm*[0:lm-1];

M=fft(m);

figure(3)

subplot(2,1,1);plot(tm,m);

grid on;

title('IIR滤波器滤波后波形');

xlabel('t');

ylabel('幅度');

subplot(2,1,2);plot(wm,abs(M));

title('IIR滤波器滤波后频谱');

xlabel('w/pi');

ylabel('幅度');

%矩形窗和布莱克曼窗

N=33;wc=0.3*pi;%基于经验的指标，其中N为理想低通滤波器阶数，wc为截止频率

hd=ideal(N,wc);%调用理想低通滤波器函数

w1=boxcar(N);%产生各种窗函数

w2=blackman(N);

h1=hd.*w1';%加窗设计各种FIR滤波器

h2=hd.*w2';

th1=(0:32)/fs;

th2=(0:32)/fs;

M=21184;

fh1=fft(h1,M);%矩形窗频谱函数

w=2/M*[0:M-1];

fh2=fft(h2,M);%布莱克曼窗频谱函数

figure(4)

subplot(2,2,1);plot(th1,h1)

title('矩形窗时域');

subplot(2,2,2);plot(w,abs(fh1));

title('矩形窗频域');

subplot(2,2,3);plot(th2,h2);

title('布莱克曼窗时域');

subplot(2,2,4);plot(w,abs(fh2));

title('布莱克曼窗频域')

%滤波处理

y6=conv(h1,y4);%用矩形窗对调制后信号进行滤波

f6=fft(y6);

n4=length(f6);

ty6=(0:n4-1)/fs;

w3=2/n4*[0:n4-1];

%sound(y6,fs);

figure(5)

subplot(2,2,1);plot(ty6,y6);

title('矩形窗滤波后音频')

subplot(2,2,2);plot(w3,abs(f6));

title('矩形窗滤波后频谱')

y7=conv(h2,y4);%用布莱克曼窗对调制后的信号进行滤波

f7=fft(y7);

n5=length(f7);

ty7=(0:n5-1)/fs;

w4=2/n5*[0:n5-1];

%sound(y7,fs);

subplot(2,2,3);plot(ty7,y7);

title('布莱克曼窗滤波后音频')

subplot(2,2,4);plot(w4,abs(f7));

title('布莱克曼窗滤波后频谱')

 实验结果与分析：

由调制信号f(t)恢复原始信号g(t)的过程称为解调。这里，cos(w0*t)信号是接收端的本地载波信号，它与发送端的载波同频同相。f(t)与cos(w0*t)相乘的结果使频谱F(W)向左、右分别移动+w0、-w0（并乘以系数1/2），得到 

g0(t)=1/2*g(t)+1/2*g(t)*cos(w0*t)

G0(w)=1/2*G(w)+1/4*[G(w-2*w0)+G(w+2*w0)]

再利用一个低通滤波器，滤除在频率为2*w0附近的分量，即可取出g(t),完成解调。

窗函数法的设计步骤：

a、给定所要求的频率响应函数Hd（w）；

b、求出hd(n)；

c、有过渡带宽及阻带最小衰减的要求，选定窗w(n)的形状以及N的大小，一般N要通过几次试探才能确定；

d、求得所设计的FIR滤波器的单位抽样响应：h(n)=hd(n)*w(n)；

e、求H(w),检验是否符合设计要求，如不满足，则需重新设计。

函数格式说明：

[n,wn]=buttord(wp,ws,rp,rs)：n为滤波器的阶数，wn为其截止频率，wp为通带截止频率，ws为阻带截止频率，rp为通带最大衰减,，rs为阻带最小衰减。

[a,b]=butter(n,rp)：a，b为所设计的滤波器传输函数分子和分母的系数。

[h,w]=freqz(a,b):求滤波器的频率响应h，w为频域坐标。

1、由调制后的信号频谱知调制过程实质为信号频域的搬移，而解调则为调制的逆过程；

 根据解调后信号的频谱，设置巴特沃斯IIR低通滤波器的通带截止频率为0.5，阻带截止 率为0.50，阻带和通带衰减分别为1db和15db;

2、由滤波器滤波后波形频谱知，与原信号频谱相比原信号频谱中的一些高频成分也被滤除，但总体来说滤波器所选指标还是能够较好地还原信号信息，通过播放滤波后音频也可以说明这一点；

3、对于窗函数的设计，其指标则是基于频谱图和过往经验的综合，为了便于比较两种窗函数滤波的效果，两种窗函数选用相同的指标；实验中设置窗函数的阶数N=33，其截止频率Wc=0.3*pi(比巴特沃斯IIR低通滤波器截止频率略大)；通过比较两窗函数的频谱和滤波后的频谱知，矩形窗通带与阻带的过渡部分比较陡，对于高频信息滤除的比较彻底；而布莱克曼窗的过渡部分则相对比较平缓，相应的其对于高频信息的滤除就没有矩形窗那么彻底（相对而言）；比较采用矩形窗和布莱克曼窗的频率特性图可以看出：最小阻带衰减只由窗形状决定，而不受阶数N的影响；而过渡带的宽度窗的形状有关；同时，布莱克曼窗的过滤带宽、旁瓣峰值和主瓣宽度均大于矩形窗的过滤带宽、旁瓣峰值和主瓣宽度；分别播放两种窗函数滤波后的声音，通过布莱克曼窗的声音比通过矩形窗的声音稍微响亮一些（不是很明显），这也印证了上面对于两种窗函数的分析。

5、音乐信号的滤波去噪

①给原始音乐信号叠加幅度为0.05，频率为3kHz、 5kHz、8kHz的三余弦混合噪声，观察噪声频谱以及加噪后音乐信号的音谱和频谱，并播放音乐，感受噪声对音乐信号的影响；

②给原始音乐信号叠加幅度为0.5的随机白噪声（可用rand语句产生），观察噪声频谱以及加噪后音乐信号的音谱和频谱，并播放音乐，感受噪声对音乐信号的影响；

③根据步骤①、②观察到的频谱，选择合适指标设计滤波器进行滤波去噪，观察去噪后信号音谱和频谱，并播放音乐，解释现象。 

源程序:

clc;clear;close;

[y,fs,bit]=wavread('I do片段');

y0=y( : ,1);

l=length(y0);

 %加三余弦混合噪声

t0=(0:l-1)/fs;

d0=[0.05*cos(2*pi*3000*t0)]';

t1=(0:l-1)/fs;

d1=[0.05*cos(2*pi*5000*t1)]';

t2=(0:l-1)/fs;

d2=[0.05*cos(2*pi*8000*t2)]';

noise=d2+d1+d0;

y1=y0+noise;

%sound(y1,fs);

a=length(noise);%绘制三余弦噪声音频图

wa=2/a*[0:a-1];

Noise=fft(noise);

figure(1)

subplot(2,2,3);

plot(noise(1:150)); 

xlabel('时间（s）')

ylabel('幅值')

title('三余弦信号音谱')

subplot(2,2,1);%绘制三余弦噪声频谱图

plot(wa,abs(Noise));

grid on; 

xlabel('W')

title('噪声频谱')

w0=2/l*[0:l-1];%绘制加噪信号音频

Y1=fft(y1);

subplot(2,2,4)

plot(w0,abs(Y1));

grid on; 

xlabel('W')

title('加噪信号频谱')

ly1=length(y1);

ty1=(0:ly1-1)/fs;

subplot(2,2,2);plot(ty1,y1);

xlabel('时间（s）')

ylabel('幅值')

title('加噪信号音谱')

m=rand(l,1)-0.5; %产生幅度为0.5的随机信号

lm=length(m);

y2=m+y0;%将噪声信号与原声音信号叠加

%sound(y2,fs);

wm=2/lm*[0:lm-1];

M=fft(m);

figure(2)

subplot(2,2,3);

plot(m(1:150))

xlabel('时间（s）')

ylabel('幅值')

title('白噪信号音谱')

subplot(2,2,1);

plot(wm,abs(M));

grid on; 

xlabel('W')

title('噪声频谱')

l=length(y2);

ty2=(0:l-1)/fs;

w=2/l*[0:l-1];

Y2=fft(y2);

subplot(2,2,4)

plot(w,abs(Y2));

grid on; 

xlabel('W')

title('加噪信号频谱')

subplot(2,2,2);plot(ty2,y2);

xlabel('时间（s）')

ylabel('加噪信号幅值')

title('加噪信号音谱'); 

%设计滤波器进行滤波去噪

[N1,wc1]=buttord(0.04,0.13,1,24);%确定低通滤波器的阶数和截止频率；

%[N1,wc1]=buttord(0.04,0.20,1,18);（针对随机白噪声的指标）

[b,a]=butter(N1,wc1);           %确定低通滤波器分子分母系数

[H,W]=freqz(b,a);

figure(3);

subplot(3,1,1);

plot(W,abs(H));%滤波器频谱

xlabel('w/pi')

ylabel('幅度k')

title('IIR滤波器滤波波形'); 

grid on;

m=filter(b,a,y1);%用滤波器滤除三余弦噪声

%sound(m,fs);

lm=length(m);%滤波后信号长度

tm=(0:lm-1)/fs;%设置横坐标

subplot(3,1,2);

plot(tm,m);%绘制滤波后的波形

xlabel('t(s)')

ylabel('信号幅值')

title('去噪后信号波形');      

k=fft(m);  %滤波后的波形做离散傅里叶变换

w=2*[0:length(k)-1]/length(k);

subplot(3,1,3)

plot(w,abs(k));

xlabel('w/pi')

ylabel('幅度k')

title('IIR滤波器滤波后信号频谱');    

实验结果与分析：

信号的去噪即根据信号噪声频谱的分布范围，从而采用不同的滤波器对其进行

滤波。常用的滤波器有IIR巴特沃斯低通滤波器，FIR窗函数滤波器。

Rand函数介绍：rand函数产生由在（0，1）之间均匀分布的随机数组组成的数组。Y = rand(n) 返回一个n x n的随机矩阵。如果n不是数量，则返回错误信息。

Y = rand(m,n) 或 Y = rand([m n]) 返回一个m x n的随机矩阵。

Y = rand(m,n,p,...) 或 Y = rand([m n p...]) 产生随机数组。

Y = rand(size(A)) 返回一个和A有相同尺寸的随机矩阵。1、rand(3)*-2 ：rand（3）是一个3*3的随机矩阵（数值范围在0~1之间）然后就是每个数乘上-2； 2、 用matlab随机产生60个1到365之间的正数  ：1+fix（365*rand（1，60））；3、用rand函数随机取100个从-1到1的数x1，x2，...：x = rand(1,100) * 2 - 1

1、通过观察三余弦噪声和随机白噪声的频谱知，三余弦噪声的频谱为三条平行于y轴的线段；而随机白噪声则为连续的不规则的波形；加三余弦噪声后的音乐片段听着带有连续的蜂鸣声，而加随机白噪声的音乐片段听起来则是那种典型的噪声，并且也是不间断的；

2、由于对巴特沃斯滤波器有更深入的了解，且掌握的比较熟练，故在此处均采用巴特沃斯滤波器进行滤波去噪。参数选择如下：

三余弦去噪：由分析其频谱可知，三余弦混合噪声信号的三条频谱线分别加在频率为w1=0.136π，w2=0.227π和w3=0.363π处，所以选取wp=0.04*pi为通带截止频率，ws=0.13pi为其阻带截止频率，rp=1db为其通带最大衰减,，rs=24db为阻带最小衰减。.

随机白噪声去噪：由于随即白噪声信号是加在了整个时域和频域内，为了滤波的更彻底，所以选取较小的通带截止频率和阻带截止频率，并选取较大的组带最小衰减。wp=0.04pi为通带截止频率，ws=0.20pi为其阻带截止频率，rp=1db为其通带最大衰减,，rs=18db为阻带最小衰减。

用对应滤波器对加三余弦噪声的音频过滤后，原音乐片段的高频信息也被滤除了一部分，播放滤波后的音乐，音乐信号中仍然夹杂着少许蜂鸣声（当然比滤波前弱了许多）；

用另一对应滤波器对加了随机白噪声音频过滤后，原音乐片段的频谱中也有相当一部分高频信息被滤除，并且播放滤波后的音乐片段时也夹杂着少许噪声；之所以无法完全滤除加在音乐片段上的噪声，一部分是因为仍有少量低频噪声叠加在音乐信号的低频部分，这是低通滤波器所无法滤出的，另一部分是因为所设计滤波器指标没有达到理想值；

6音乐信号的幅频滤波及相频分析

① 设计低通滤波器，滤除原始音乐信号的高频信息，观察滤波前后的幅度频谱， 并比较滤波前后音乐效果，感受高频信息对音乐信号的影响；

②设计高通滤波器，滤除原始音乐信号的低频信息，观察滤波前后的幅度频谱， 并比较滤波前后音乐效果，感受低频信息对音乐信号的影响；

③选取两段不同的音乐信号，分别将其幅度谱和相位谱交叉组合构成新的音乐信号，播放并比较组合后音乐与原始音乐，感受相频信息对音乐信号的影响。

源程序：6-1、2

clear all;clc

[y,fs,bit]=wavread('I do片段');

size(y)%查看读取信号的声道类型

y1=y(: ,1);%对信号进行分列处理

n=length(y1);%求信号y1的的长度

t1=(0:n-1)/fs;

f1=fft(y1);%对y1进行fft谱分析

w=2/n*[0:n-1];%w为连续频谱的数字角频率横坐标

%sound(y,fs);%播放音乐信号

figure(1)

subplot(2,1,1);plot(t1,y1);

title('音乐信号的波形');

xlabel('t');

ylabel('y1');

subplot(2,1,2);plot(w,abs(f1));

title('音乐信号的频谱');

xlabel('w');

ylabel('f1');

%用IIR滤波器滤波（低）

[n2,wc2]=buttord(0.15,0.20,1,15);%确定低通滤波器的阶数和截止频率；

[B2,A2]=butter(n2,wc2); %确定低通滤波器分子分母系数

[H,W]=freqz(B2,A2);

figure(2)

subplot(3,1,1);

plot(W,abs(H));%低通滤波器波形

xlabel('w/pi')

ylabel('幅度k')

title('滤波器（低通）频谱'); 

grid on;

m2=filter(B2,A2,y1);%滤波

lm2=length(m2);

tm2=(0:lm2-1)/fs;

subplot(3,1,2)

plot(tm2,m2);

xlabel('n')

ylabel('信号幅值')

title('低通滤波器滤波后信号波形');   

%sound(m2,fs);

k2=fft(m2);                        %滤波后的波形做离散傅里叶变换

l2=length(k2);                        

w2=2*[0:l2-1]/l2;

subplot(3,1,3);

plot(w2,abs(k2));

xlabel('数字角频率w')

ylabel('幅度')

title('低通滤波器滤波后信号频谱');  %解调滤波后的频谱

%用IIR滤波器滤波（高）

[N,WC]=buttord(0.15,0.20,1,15);%确定高通滤波器的阶数和截止频率；

[B,A]=butter(N,WC,'high');   %确定高通滤波器分子分母系数

[H1,W1]=freqz(B,A);

figure(3)

subplot(3,1,1);

plot(W1,abs(H1));%高通滤波器频谱

xlabel('w/pi')

ylabel('幅度k')

title('滤波器（高通）频谱'); 

grid on;

m=filter(B,A,y1); %滤波

lm=length(m);

tm=(0:lm-1)/fs;

subplot(3,1,2);

plot(tm,m);

xlabel('n')

ylabel('信号幅值m')

title('高通滤波器滤波后信号波形');      

%sound(m,fs);

k=fft(m); %滤波后的波形做离散傅里叶变换

l2=length(k);                        

w2=2*[0:l2-1]/l2;

subplot(3,1,3);

plot(w2,abs(k));

xlabel('数字角频率w')

ylabel('幅度k')

title('IIR高通滤波器滤波后信号频谱');  

实验结果与分析：

1、用低通滤波器滤波后，声音变小了，但听起来很刺耳，失真很严重，无法辨别歌声；

2、用高通滤波器滤波后，声音听起来很低沉，听起来就好像有杂音一样，但基本还能辨别原歌曲音调；

6-3选取两段音乐信号，分别将其幅度谱和相位谱交叉组合构成新的音乐信号

clear all;clc

[a,fs1,bit1]=wavread('I do片段');

[b,fs2,bit2]=wavread('风声片段');

size(b)%查看读取信号的声道类型

y2=b(: ,1);%对信号进行分列处理

n2=length(y2);%求信号y2的的长度

t2=(0:n2-1)/fs2;

f2=fft(y2);

w2=2/n2*[0:n2-1];

%wavplay(y2,fs2);

size(a)%查看读取信号的声道类型

y1=a(: ,1);%对信号进行分列处理

n1=length(y1);

t1=(0:n1-1)/fs1;

f1=fft(y1);

w1=2/n1*[0:n1-1];%w为连续频谱的数字角频率横坐标

%wavplay(y1,fs1);

Fy1=abs(f1);%音乐1的幅度

Ay1=angle(f1);%音乐1的相位

Fy2=abs(f2);%音乐2的幅度

Ay2=angle(f2);%音乐2相位

F1=Fy1.*exp(j*Ay2);%音乐1的幅度与音乐2的相位交叉组合

X1=ifft(F1);

n3=length(X1);

tx1=(0:(n3-1))/fs1;

w3=2/n3*[0:n3-1];

%wavplay(real(X1),fs1);%播放交叉组合后的音乐

F2=Fy2.*exp(j*Ay1);%音乐2的幅度与音乐1的相位交叉组合

X2=ifft(F2);%音乐2的幅度与音乐1的相位交叉组合后的信号函数

n4=length(X2);

tx2=(0:(n4-1))/fs2;

w4=2/n4*[0:n4-1]

%wavplay(real(X2),fs2);%播放交叉组合后的音乐

figure(1)

subplot(2,1,1);plot(t1,y1);%绘制音乐信号1的波形

title('音乐信号1的波形');

xlabel('t');

ylabel('y1');

subplot(2,1,2);

plot(w1,abs(f1));%绘制音乐信号1的频谱

title('%音乐信号1的频谱');

xlabel('w');

ylabel('f1');

figure(2)

subplot(2,1,1);plot(t2,y2);%绘制音乐信号2的波形

title('音乐信号2的波形');

xlabel('t');

ylabel('y2');

subplot(2,1,2);plot(w2,abs(f2));

title('音乐信号2的频谱');

xlabel('w');

ylabel('f2');

figure(3)

subplot(2,1,1);plot(tx1,X1);

title('音乐1的幅度与音乐2的相位交叉组合后的波形');

xlabel('t');

ylabel('X1');

subplot(2,1,2);plot(w3,abs(F1));

title('音乐1的幅度与音乐2的相位交叉组合后的频谱');

xlabel('w');

ylabel('F1');

figure(4)

subplot(2,1,1);plot(tx2,X2);

title('音乐2的幅度与音乐1的相位交叉组合后的波形');

xlabel('t');

ylabel('X2');

subplot(2,1,2);plot(w4,abs(F2));

title('音乐2的幅度与音乐1的相位交叉组合后的频谱');

xlabel('w');

ylabel('F2');

实验结果与分析：

1、音乐1的幅度和音乐2的相频交叉组合，播放产生的新的音乐信号，还能听得出音乐2的音调，并且与原音乐2的音调相比夹杂了一些噪音；

2、音乐2的幅度和音乐1的相频交叉组合，播放产生的新的音乐信号，还能听得出音乐1的音调，并且与原音乐1的音调相比夹杂了一些噪音；

3、音乐信号的相频决定了不同频率成分震动的初相位，而幅频决定了不同频率成分震动的幅度大小，综合1,2知相频在音乐信号起主要作用（就音乐信号的旋律而言）。

四、学习matlab的心得体会：

通过近两个月的学习，我觉以下几点对于学好matlab很重要：兴趣、悟性和坚持。

 1．多动手写程序、调试，先有量变后有质变

下面是一些常见错误提示信息

1.Subscript indices must either be real positive integers orlogicals

中文解释：下标索引必须是正整数类型或者逻辑类型

出错原因：在访问矩阵（包括向量、二维矩阵、数组，下同）的过程中，下标索引要么从 0 开始，要么出现了负数。注：matlab的语法规定矩阵的索引从 1 开始，这与 C 等编程语言的习惯不一样。

解决办法：自己调试一下程序，把下标为 0 或者负数的地方修正。

2.Undefined function or variable "a"

中文解释：函数或变量 a 没有定义

3.Input argument "x" is undefined

中文解释：输入变量 x 没有定义

4.Matrix dimensions must agree

Inner matrix dimensions must agree

中文解释：矩阵的维数必须一致

出错原因：这是由于运算符(= + - / * 等)两边的运算对象维数不匹配造成的，典型的出错原因是错用了矩阵运算符。matlab通过“.”来区分矩阵运算和元素运算

5.Function definitions are not permitted at the prompt or inscripts

中文解释：不能在命令窗口或者脚本文件中定义函数

出错原因：一旦在命令窗口写 function c = myPlus(a,b)，此错误就会出现，因为函数只能定义在 m 文件中

6. 1） X must have one or two columns

   2）Vectors must be the samelengths

中文解释：

1. X 必须是 1 或者 2 列

2. 向量长度必须一致

7.One or more output arguments not assigned during call to'...'

中文解释：在调用...函数过程中，一个或多个输出变量没有被赋值

8. Error using ==> mpower

Matrix must be square

中文解释：错误使用mpwoer函数，要求矩阵必须是方阵

9.Explicit integral could not be found.

中文解释：显式解没有找到

10.Index exceeds matrix dimensions.

Attempted to access b(3,2); index out of bounds becausesize(b)=.

中文解释：索引超出矩阵的范围

11.In an assignment   A(I) =B, the number of elements in B and I must be the same

中文解释：在赋值语句 A(I) = B 中，B 和 I 的元素个数必须相同

12.To RESHAPE the number of elements must not change

中文解释：矩阵变换时，变换前和变换后的总元素不能改变

正确的理解了错误提示信息的含义，程序的修改就会变得事半功倍。

2．善于利用MATLAB的帮助

一遇到什么问题,通常我的第一反应是:help，就先说说自己对help的一些常用方法吧。

1）命令窗口直接敲“help”，你就可以得到本地机器上matlab的基本的帮助信息。

2）对于某些不是很明确的命令，只知道大体所属范围，譬如说某个工具箱，直接在命令窗口中敲入

Help toolboxname，就可以得到本工具箱有关的信息：版本号，函数名等。

3）知道函数名，直接用help funname就可以得到相应的帮助信息。

3．善于向别人学习，多看别人写的代码并消化

三人行，必有我师，一个人很难什么都精通，取长补短是最快的进步方法。

五、探究与创新：

1、IIR数字滤波器和FIR数字滤波器的比较

答：IIR滤波器：相位一般是非线性的，不一定稳定，不能使用fft做快速卷积，一定是递归结构，对频率分量的选择性好（零极点可同时起作用），相同性能下阶次较低，有噪声反馈，噪声大，运算误差大，有可能出现极限环振荡，设计时有大量图表可查，方便简单，主要用于设计分段常数的标准低通，带通，高通，带阻和全通滤波器。

     FIR滤波器：相位可以做到严格线性，一定是稳定的，信号通过系统可采用快速卷积，主要是非递归结构，也可含递归结构，选择性差，相同性能下阶次高，噪声小，运算误差小，不会出现极限环振荡，可设计正交变换器、微分器、线性预测器、线性调频器等各种网络，适用范围广。

     其中，FIR滤波器的最大好处是稳定、线性相位和广泛的适用范围，而它的最大

缺点是阶数高，从而带来时延大、存储单元多等问题。例如用频率抽样法设计阻带衰耗为－20dB的FIR DF需33阶，用双线性法设计同样指标的切比雪夫IIR DF仅需4～5阶。因此，在一些对时延有严格的场合就不得不考虑用IIR滤波器。语音信号对相位的非线性不很敏感。数据和图象信号则往往对滤波器提出线性相位的要求，这就是为什么FIR用得越来越广的原因。总之，IIR和FIR各有特点，在应用时要根据各方面的指标，综合考虑加以选择。

2、音乐信号的音调与信号的什么特征有关？

答：音调主要由声音的频率决定。对一定强度的纯音，音调随频率的升降而升降；对一定频率的纯音、低频纯音的音调随响度增加而下降，高频纯音的音调却随响度增加而上升。音调的高低还与发声体的结构有关，因为发声体的结构影响了声音的频率。所以音乐信号的音调与信号的频率有关。音调还与声音持续的时间长短有关。非常短促(毫秒量级或更短)的纯音，只能听到像打击或弹指那样的“喀嚓”一响，感觉不出音调。持续时间从10 毫秒增加到50 毫秒，听起来觉得音调是由低到高连续变化的。超过50 毫秒，音调就稳定不变了。

3、音乐信号的音色与信号的什么特征有关？

答：音乐信号的音色是指对声音音质的感觉特性，也是一种声音区别于另一种声音的特征品质。音色的不同取决于不同的泛音，每一种乐器、不同的人以及所有能发 

声的物体发出的声音，除了一个基音还有许多不同频率的泛音伴随，正是这些泛音决定了其不同的音色，使人能辨别出是不同的乐器甚至不同的人发出的声音。而泛音的不同指的是频谱的不同，所以音乐信号的音色与信号的频谱或频率特性有关。

4、两种不同音色的音乐信号叠加混叠后，为何人耳还可以分辨？

答：人耳听觉系统实际上就是一个音频信号处理器，可以完成对声音信号的传输、转换以及综合处理的能力。而音乐信号的音色与音乐信号的频谱分布有关，不同音色的音乐信号的频谱当然不会相同，当混叠后的音乐信号进入人耳的听觉系统后，人耳通过对混叠信号的频谱处理，可以将具有不同频谱的声音信号分离开，也就达到了分辨不同音色的音乐信号的目的。

5、音乐信号的幅度与相位特征对信号有哪些影响？

答：音乐信号的相位发生改变会引起信号时域特性发生改变；音乐信号的幅度发生变化会影响信号的频谱幅度。