当前位置:首页 >> 建筑/土木 >>

[一] 课程设计目的 在PC机平台上用MATLAB语言设计


三、音乐分析系统

[一] 课程设计目的 在 PC 机平台上用 MATLAB 语言设计一个音乐分析系统,对一段音乐做 DFT,通过频谱分析 的方法来确定音符及其节拍。 [二] 课程设计的内容和要求 一、窗函数分析 对一段信号做频谱分析,必须进行加窗处理。因此,窗函数类型的选择非常重要。用于频 谱分析和滤波器设计中的窗函数有许多种类型,本题目首先分析不同类型的窗

函数,对它们的 频率响应进行评价。 设 x(n ) 为一长序列, w(n ) 为长度为 N 的窗函数, n = 0,1, 得

, N ? 1 ,用 w(n ) 乘以 x(n ) ,


x N (n ) = x(n )w(n )

x N 为 N 点序列,其离散傅立叶变换为

X N e jω =

( )

1 2π

π ω (ω θ ) ∫ π X (e )W (e )dθ
j j ? ?



由此可以看出, 窗函数 w(n ) 不仅影响原信号在时域的形状, 也影响了其在频域的形状。 w(n ) 若 为矩形窗,即 w(n ) = 1, n = 0,1,

, N ? 1 ,其频谱为
N ?1 n =0

W e jω = ∑ w(n )e ? jωn =e
? jω N ?1 2

( )

? ωN ? ?ω ? sin ? ? /sin ? ? ? 2 ? ?2?



? ωN ? sin ? ? ? 2 ? 是其幅度响应。 其中 Wr (ω ) = ?ω ? sin ? ? ?2?
归一化的幅度响应如图 3.1 所示。 在图中, ≤ ω

2π 2π 的部分称为窗函数的主瓣, ω < 在 N N

后的部分称为窗函数的旁瓣。根据式⑵及 DFT 的性质可知,窗函数的主瓣宽度决定了被截短后 所得序列的频域分辨率,而旁瓣峰值有可能淹没信号频谱分量中较小的成分。因此,我们对窗 函数的总的要求是,希望它的 频谱的主瓣宽度尽可能地窄,而旁瓣的峰值尽可能地小,从而 使得信号的能量尽量集中在主瓣内。
12

0

-10

-20

度 幅

-30

-40

-50

-60

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

归 一化 频 率 , fs=1

图 3.1

矩形窗函数的归一化幅度响应(N=8)

常见的窗函数有以下几种: 1.矩形窗:

?1, 0 ≤ n ≤ N ? 1 ? w(n ) = ? 其它 ?0, ?
2.三角窗:



N ?1 ? 2n 0≤n≤ ? N ?1, 2 ? 2n N ?1 ? , w(n ) = ?2 ? ≤ n ≤ N ?1 2 N ?1 ? ?0, 其它 ? ?
3.汉宁窗:



? ? ? 2πn ?? ??, ?0.5?1 ? cos? ? ? ? N ? 1 ?? w(n ) = ? ?0, ? ?
4.汉明窗:

0 ≤ n ≤ N ?1


其它

? ? 2πn ? ?0.54 ? 0.46 cos? N ? 1 ?, 0 ≤ n ≤ N ? 1 ? ? w(n ) = ? ?0, 其它 ?
5.布兰克曼窗:



13

? ? 2πn ? ? 4πn ? ?0.42 ? 0.5 cos? N ? 1 ? + 0.08 cos? N ? 1 ?, ? ? ? ? w(n ) = ? ?0, ?
6.凯泽窗:
2 ? 2n ? ? ? I 0 ? β 1 ? ?1 ? ? ? ? ? N ? 1? ? ?, w(n ) = ? I 0 (β )

0 ≤ n ≤ N ?1


其它

0 ≤ n ≤ N ?1



? 其中 I 0 [] 是修正零阶贝塞尔函数,而 β 是一个与 N 有关的参数,通过选择不同的 β 值,
可以得到不同的过渡带宽度和近似最优的阻带衰减。 这一部分的内容和要求如下: ㈠. 窗函数主瓣宽度的分析与观测。 测量不同类型的窗函数的频率响应的主瓣宽度, 并与矩 形窗做一比较。 用补零 FFT 来计算窗函数的离散傅立叶变换时, 为什么 FFT 的长度至少 是窗长度的 4 至 5 倍? ㈡. 窗函数旁瓣高度的分析与观测。测量三角窗、汉宁窗、汉明窗、布兰克曼窗等窗函数的 最大旁瓣高度,旁瓣高度与窗函数的长度有关吗? ㈢. 凯泽窗中的参数 β 是用来改变旁瓣高度和主瓣宽度的,请绘出不同 β 值情况下的凯泽

窗的时域波形及其幅度响应,并验证其最大旁瓣高度正比于

sin c(β )

β



㈣. 观察包含两个频率分别为 0.20 弧度和 0.25 弧度的单频信号的混合信号的频谱, 为了清 楚地显示出频谱中的这两个谱峰,对于不同的窗函数,分析各需要多少样本点。 要求使用 MATLAB 语言来分析窗函数,为每一部分设计内容编写出相应的 MATLAB 程序(.m 文件) ,绘出对应的波形图,并做出详细地分析。 附 1. 窗函数分析的 MATLAB 参考程序代码 N=8; 窗函数的长度 w=boxcar(N); 这是一个产生矩形窗的内建函数 W=fft(w,N*128); 计算长度为N*128的补零FFT WdB=db(W,80); 对幅度响应进行归一化处理,同时,将最大衰减统一取为-80dB plot((0:1:N*128-1)/(N*128-1),WdB); 绘出波形图 xlabel('归一化频率, fs=1');ylabel('幅度'); 二、音乐分析 这一部分的任务是,对一段音乐信号进行分析,求出其频率响应,用等高线或灰度图的形 式画出其频谱图,利用峰值检测算法找出频率中的峰值,以此来判断每个音符的持续时间,确 定音乐的节拍。 主要内容有:

14

㈠.音乐文件的读取与播放 在本课题设计中,采用.wav 文件来做分析和测试。在本设计中,大多数音乐文件的采样 频率是 f s = 11025Hz 。在 MATLAB 语言中,音乐文件的读取和播放命令是: [s,fs]=wavread('cscal_np.wav'); 读波形文件cscal_np.wav,并获取其采样频率 sound(s,fs); 播放波形文件 musicg_x; 启动音乐GUI(图形用户接口),以便可观察到音符 其中,函数 musicg_x()及其所有有关的 M 文件和 MAT 文件已在课题设计包中提供给大家,不 须要重新编写。在运行 musicg_x()后,输入音调和音符周期,即可观察到五线谱和声谱图, 同时可以播放所谱写好的歌曲。 ㈡.计算音乐信号的频谱图 一段音乐是持续时间很长的信号,因此,音乐信号的频谱必须分段计算,即,采用滑动窗 DFT 方法来计算音乐信号的频谱。滑动窗 DFT 法又称短时傅立叶变换 STFT,这是时间有关的傅 立叶变换,即,沿着时间轴滑动窗,每一次只取信号的一部分。序列 x(n ) 的 STFT 定义为:

X STFT e jω , n =
其中, w(n ) 是窗函数序列。

(

) ∑ x[n ? m]w[m]e
∞ m = ?∞

? jωm



设 DFT 的长度为 N ,这是计算分辨率;而 R 是窗函数的长度,这是物理分辨率。其中, N 和 R 的关系是: N ≥ R 。则,用 DFT 来计算 STFT,得

X STFT (k , n) = X STFT e , n |

(



)

ω=

2πk N

= ∑ x[n ? m]w[m]e
m =0

R ?1

?j

2πkm N



窗函数的功能是截取信号的一部分用于分析,并确保信号 x(n ) 所抽取的部分是近似平稳 的。由于音乐信号是一个非平稳信号,因此,窗函数的长度 R 的选择必须小心处理。一般来说, 窗越宽,参与运算的时间样本越多,STFT 的频率分辨率就越好。另一方面,窗越短,STFT 的 时间分辨率就越好,因为此时所截取的信号其平稳性越好。因此,一个较短的窗可提供一个宽 带频谱图,而一个较长的窗则得到一个窄带频谱图。 通常,窗与窗之间要有重叠,这样才能确保所有输入信号都能参与分析和处理。设 L 为窗 与窗之间的距离,则为了正确地分析音乐信号频谱,必须使得

L≤R≤N
用 STFT 计算音乐信号的频谱图的例子程序如下:
[s,fs]=wavread('C_major.wav'); nwin = input('窗函数的长度 ='); noverlap = input('窗中心间距 = '); nfft = input('FFT的长度 = '); window=hamming(nwin); % 读取波形文件s及其采样速率fs % 窗函数的长度 % 相邻数据段之间的中心间距 % FFT的长度



% 选择汉明窗,也可以是其它任意窗函数

15

specgram(s,nfft,fs,window,noverlap);

% 计算频谱图

其中,C_major.wav是一段C大调音乐,其五线谱如下图所示:

图 3.2

音乐信号 C_major.wav 的五线谱

所对应的频谱图则如图3.3所示,其中,窗宽R=1024,窗重合区间L=256,FFT长度N=1024。
2000 1800 1600 1400 1200 y c n e u q e r F 1000 800 600 400 200 0

0

1

2

3

4 Time

5

6

7

图 3.3 音乐信号 C_major.wav 的频谱图 ㈢.谱峰检测 音乐信号的节拍可通过频域找峰的方法来确定。首先计算出音乐信号的频谱,然后编写出 一个峰值检测算法程序,在音乐信号的频域中找峰,以此来判断每个音符的持续时间,从而确 定音乐的节拍。 从图3可以看出,这段C大调音乐的每个音符都是单音调信号,没有谐波成分。因此,可以 非常容易地从频谱图中找出每个音符的起始和结束时间。如果有谐波成分存在,则需要先对谐 波分量进行判断,然后才能确定每个音符的持续时间。以C调音乐cscal_np.wav为例,其频谱 图如图3.4所示,其中,窗宽R=1024,窗重合区间L=256,FFT长度N=1024。

16

5000

4000

y c n e u q e r F

3000

2000

1000

0

0

0.5

1

1.5

2

2.5 Time

3

3.5

4

4.5

5

图 3.4

音乐信号 cscal_np.wav 的频谱图

从图 4 可以看出,该音乐信号的谐波分量是非常强的,因此,在找峰过程中,要仔细区分 主频分量和谐波分量。 频率分辨率取决于窗的宽度,窗越宽,频率分辨率越高,相应的时间分辨率越低;相反, 窗越窄,频率分辨率越低,而时间分辨率就越高。我们把利用较宽的窗得到的频谱图称为窄带 频谱图,而把利用较短的窗得到的频谱图称为宽带频谱图。 这一部分的要求如下: 1. 利用音乐声谱实验程序 musicg_x(),输入一段音乐曲子(如,《两只老虎》等),观 察并分析其五线谱图和声谱图; 2. 自己编写一个音乐分析程序,该程序需要实现以下功能: ①、播放音乐信号; ②、计算频谱图(窗宽、窗重合区间、FFT 长度均可选) ; ③、利用频域峰值检测算法确定每个音符的起始和结束时间。 3. 利用自己编写的音乐分析程序,分析宽带频谱图和窄带频谱图对音乐声谱分析的影响。 并分析确定对于不同节拍的音乐信号,窗函数宽度、窗重合区间以及 FFT 长度应该怎样 选择才能使的声谱分析效果最好(时间分辨率和频率分辨率的折中) ; 4. 独奏音乐同时只有一个音符播放,而合奏音乐同时有多个音符播放。观察并分析这两种 情况下的音乐分析程序的不同结果。 附 2:峰值检测算法参考程序 % 频域找峰算法 function [peaksok]=peakfind(filename); N=8192;
% 读入文件中的一段数据

17

eval(['[f,fs]=wavread(''',filename,''',[1,2048]);']); % 对信号加窗,计算其FFT,再将幅度转换为分贝(dB)表示 F=fft(f.*kaiser(2048,10),N); Fdb=20*log(abs(F)); n=(1:N/2)*fs/N; % plot results on different scales figure(1);clf; plot(n,Fdb(1:length(n))); xlabel('频率'); ylabel('幅度(dB)'); title('原始信号频谱'); % 根据最大值的百分比来找出最大值,并且信号频谱进行剪辑,只保留幅值过半的频点 maxF=max(Fdb), Fdb=Fdb.*(Fdb>.5*maxF); %plot results on different scales figure(2);clf; plot(n,Fdb(1:length(n))); xlabel('频率'); ylabel('幅度(dB)'); title('幅值过半的频点'); figure(3);clf; plot(n(1:length(n)/2),Fdb(1:length(n)/2)); axis([1,length(n)/2,floor(.5*maxF),ceil(maxF)]); xlabel('频率'); ylabel('幅度(dB)'); % find peaks using the derivatives fd=diff(Fdb(1:N/2)); fdplus=(fd>0)-(fd<0); [ind]=find(abs(diff(fdplus))>1); peaksok=[ind*fs/N,Fdb(ind)], % 对峰值排序,寻找谐波分量 [sortf,indf]=sort(Fdb(ind)); sortf=flipud(sortf); indf=flipud(indf); sortf,indf, % search top few signals for harmonics, if at a multiple % and 10dB down then call it a harmonic harm = ones(size(indf)); for i=1:4

18

index=indf(i);freq=peaksok(index,1);i' if (freq>0) level=peaksok(index,2); level,freq ff=rem(peaksok(:,1),(freq*ones(length(indf),1))); harm=find(ff<.18|(freq-ff)<10);ff',harm', for j=1:length(harm) if (peaksok(harm(j),2)+10<level) peaksok(harm(j),1)=0; end; end; peaksok, end; end;

参考资料:
1. 《精通 MATLAB 6》 ,张航 黄攀译,清华大学出版社,2002 年 6 月 2. 《Digital Signal Processing: A Computer-Based Approach》 ,Second Edition,清华大 学出版社 2001 年影印版

19


相关文章:
MATLAB课程设计任务书
巩固并加深对 MATLAB 语言程序设计知识理解; 2、...课程设计报告要求: 1.封面、目录; 2.课程设计目的...检查学生进度, 学生不得以自己有私人电脑为借口而不...
matlab课程设计模板
课程设计辅导资料: “Matlab 语言基础及使用入门”“数字信号处理原 、 理与...几乎所有 PC 机和大型计算机 上运行,适用于 Windows、UNIX 等多种系统平台。...
matlab课程设计Word文档
大学课程设计报告 课程设计任务书一.本课程设计的目的 (1)熟练掌握 MATLAB 的...《MATLAB 语言及系统仿真课程设计》 暨南大学课程设计报告 误差分析 工具应用 ...
MATLAB课程设计
数字信号处理 课程设计报告 题目: MATLAB 余弦离散傅里叶变换 学院: 培黎工程...一. 设计目的 1.要求学生会用 MATLAB 语言进行编程,绘出所求波形,并且运用 ...
matlab课程设计 (20)
在数字信号处理平台上 (PC 机MATLAB 仿真软件系统) 进行软件仿真 设计,并进行调试和数据分析。 2. 利用 MATLAB 仿真软件进行绘图。 课程设计的目的: 1.理论...
MATLAB课程设计
一、课程设计目的 1.掌握 MATLAB 基本操作; 2.学会利用 MATLAB 进行基本...提供了一种全面解决方案,并在很大程度上摆脱了传统 非交互式程序设计语言(如...
MATLAB课程设计
MATLAB 课程设计报告题目:利用 MATLAB 实现信号幅度调制与解调 班姓学姓学 级...离散傅叶变换的一种快速算法,简称 FFT。 为了节省电脑的计算时间,实现数字信...
matlab实验设计
一.课程设计目的 1:熟悉课程设计的基本流程。 2:掌握 MATLAB 语法结构及调试方法。 3:熟悉 MATLAB 函数调用,熟练二维画图。 4:掌握 MATLAB 语言在控制方面运用...
心电信号处理课程设计-matlab
本课题的目的设计课题主要研究数字心电信号的初步...采用 Matlab 语言设计,要求分别采用两种方式进行仿真,...三、主要设备和软件 (1)PC 机一台。 (2) ...
MATLAB实验设计报告_图文
课程中较多的知识点;也可以自行学习课本上未涉及的 Matlab 程序设计线 语言高级...基于 MATLAB 的图像处理的课程设计 一、设计目的: 设计目的:综合运用 MATLAB ...
更多相关标签:
matlab课程设计目的 | c语言课程设计目的 | matlab实验目的 | matlab曲线拟合的目的 | 课程设计的目的 | 课程设计的目的和意义 | 单片机课程设计目的 | java课程设计目的 |