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

[一] 课程设计目的 在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 ,其频谱为

W e jω = ∑ w(n )e ? jωn
n =0

( )

N ?1

=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



? 是修正零阶贝塞尔函数,而 β 是一个与 N 有关的参数,通过选择不同的 β 值, 其中 I 0 []
可以得到不同的过渡带宽度和近似最优的阻带衰减。 这一部分的内容和要求如下: ㈠. 窗函数主瓣宽度的分析与观测。 测量不同类型的窗函数的频率响应的主瓣宽度, 并与矩 形窗做一比较。 用补零 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


相关文章:
数字信号处理课程设计报告1
13 1.课程设计目的及意义 (1)学会 MATLAB使用...计算机,也可以将所需要的运算编成程序,让通用计算机...FORTRAN 等语言完相 同的事情简捷得多,工具包又...
Matlab编程与系统仿真:基带传输课程设计资料_图文
Matlab编程与系统仿真:基带传输课程设计资料_计算机软件及应用_IT/计算机_专业资料...(噪声为加性高斯白噪声) 二、设计目的 1.综合应用《Matlab 编程与系统仿真》 ...
实验一:用matlab的linprog求解简单线性规划问题
实验目的: (1)使学生在程序设计方面得到初步的训练...学习 Matlab (C 或 VB)语言进行程序设计中一些 ...二、实验用仪器设备、器材或软件环境 计算机, Matlab...
利用MATLAB仿真软件系统结合双线性变换法设计一个数字...
摘要 Matlab 是一个矩阵设计平台,传统数字滤波器设计需要大量的计算,但是利用 Matlab 可以快速实现滤波器的设计与仿真,而且频谱分析功能强大,在数字信 号处理中发挥...
噪声产生器的MATLAB实现及性能分析——噪声带宽为1.7MH...
课程设计任务书计算机与通信工程 学院 通信工程 专业...1.7MHz 主要内容: 本课程设计目的主要是仿真噪声...1.3 设计平台设计开发平台MATLAB 中的 ...
在matlab上的的QPSK调制与解调仿真
课程设计中,系统开发平台为 Windows 2000,程序运行...11 利用 MATLAB 实现 QPSK 调制及解调 1.前言 1....Matlab 的软件仿真,只需 PC 机上安装 MATLAB 6.0 ...
MATLAB通信仿真
一、课程设计目的和要求 1.综合运用信号与线性系统,通信原理,以及 matlab 语言...硬件:PC机 2. 软件:MATLAB7.1,Windows7操作系统 三:课程设计原理及背景知识...
基于MATLAB的图像处理的课程设计(车牌识别系统)
MATLAB 既是一种直观、高效的计算机语言,同时 又是一个科学计算平台。它为数据...关键词:MATLAB,数字图像处理,车牌识别系统 2 一、课程设计目的 Matlab 技术课程...
一个matlab课程设计报告例子
界面友好, 最常用的方法就是使用图形界面,而 Matlab 是一门面向对象的 语言。...任务描述 1.课程设计目标 用 GUIDE 编写一个简易计算器和一个能自动画图的程序...
信号与系统课程设计1
目录一、 设计目的及意义...信号与系统课程设计 课程名称: 题目名称: 学院: ...此外, 利用 MATLAB 中的函数, 可以便捷的得到给定系统的冲激响应图像。 因此,...
更多相关标签: