小波变换原理与基本案例分析

科研训练报告

一：科研训练目标及意义

在深入理解经典Fourier 变换的基础上，通过系统介绍加窗Fourier 变换（Windowed Fourier Transfrom, WFT ）,小波变换（Wavelet trsnaform,WT ）等定义，理解时间分辨率，频率分辨率及，时频局部化性质等基本概念；通过编程实现基本的时频分析方法，理解时变信号中时频结构信息。

二：科研训练任务

1.理解加窗Fourier 变换的意义

2.编程实现加窗Fourier 分析几种合成信号

3.利用加窗Fourier(小波变换)分析几种典型的合成信号

4.观察分析窗函数，窗长变化，频率，分辨率的影响

5.可尝试利用加窗Fourier （小波变换）正变换分析实测语音等信号

三:报告结构

四：加窗Fourier 分析 v.1.Fourier 正变换

v.1.a 离散Fourier 正变换数学描述

有限长序列通过离散傅里叶变换，并通过快速傅里叶变换算法FFT 即得：

X (k )=∑

x (j )W ?

(j?1)(R?1)

?

j=1

其中：W ?=e ?2?π?j/N v.1.b.Fourier 正变换Matlab 实现 被分析信号为N 维列向量；

x =[x (0),x (1),x (2),x (3)?x (n ?1)]T

变换因子为：

[W 0?W 0

???

W 0

?W

(N?1)(N?1)

] DFT 结果为

X =[X (0),X (1),X (2),X (3)?X (n ?1)]T

以矩阵形式表示DFT ： X N =W N x N （具体代码参见代码code erst ） v.1.c.Fourier 正变案例分析总结 step1 检测代码的正确性

当输入一个离散的冲激函数的Delta 函数，其DFT 变换的幅度谱为一等幅，幅度为1的白色谱。（参见图a.3.1） Step2 分析信号特性

当输入一个正弦函数时，运行结果如下：幅度谱为一对称的离散冲激信号，冲激的横坐标为其频率，纵坐标为其幅度。由观测可看出在主

频的附近有明显的频谱泄露。（参见图a.3.2）

当输入等幅的正弦波叠加则运行结果显示，两个正弦波的幅值并不相同，存在明显的栅栏现象。（参见图a.3.3） 当输入一分时频信号时，该信号表达式为：

y ={

t <0.5, (sin (2?pi ?50?t1))t ≥0.5, (sin (2?pi ?120?t1))

其快速傅里叶变换幅频特性显示如图a.3.4，与图a.3.3相比，其信号频率仍然为50，120hertz 。并无变化。

图v.a.3.1 输入信号为序列长度为32的离散冲激信号

Xot(1)与Xot(2)的区别在于其相位不同，后面是经过DFT 变换的幅度谱和相位谱

横坐标为点数，纵坐标为幅值和相位。

图v.1.c.1 输入信号为sin(2*pi*50*t)，横坐标为模拟频率，纵坐标为其幅度值。

图v.1.c.2 输入信号为sin(2*pi*50*t)+sin(2*pi*120*t) 横坐标为模拟频率，纵坐标为其幅度值。

图v.1.c.3输入信号为一分时频信号，在t<0.5时，其表达式为y1=sin(2*pi*50*t1) t>0.5时，其表达式sin(2*pi*120*t2)。横坐标为模拟频率，纵坐标为幅值。 横坐标为模拟频率，纵坐标为其幅度值。

2.加窗Fourier 分析

2.a 加窗Fourier 正变换数学描述

通过参见图a.3.1和图a.3.3和图a.3.4对比分析可见，傅里叶变换是一种针对-(?∞,∞)的全局变换。

而加窗Fourier 变换的基本思想是在时域内进行分段，用位置不同的窗函数与原信号相乘。在时域内，时窗函数一般选取具有能量局部化的函数，选定一个基本窗函数后，通过将窗函数中心沿时间轴平移，得到一组窗函数，使得窗函数很好的覆盖时间轴。因此，时间窗等效提取了源信号的不同时间内的信号而时间带外的信号衰减为零。

加窗Fourier 变换认为在窗内信号是准平稳的，可采用平稳信号的分析方法，如频谱分析。

最早的Gabor 的定义如下;

设，g (t )∈L 2(R )即0<∫|g (t )|2?t ∞

?∞<∞，且g (t )为实对称函数，则信号的窗口傅里叶变换定义为：

G f (ω,b)=∫f (t )g (t ?b )e ?jωt ?t ∞

?∞

* 参数解释：g (t )为基本床函数，g (t ?b )为被平移到t=b 处的窗函数， G f (

)

=∫f ?e ?jωt ?t ∞

?∞

，即为古典傅里叶变换。 由其数学描述可见G f (ω,b)分析了t=b 附近的信号的信息。 2.b. 加窗Fourier 正变换Matlab 实现

通过自定义函数实现加窗傅里叶变换，函数的重要参数为测试信号x,窗口长度Nw ，窗函数类型。

一般选取分时频信号作为窗口傅里叶变换的测试信号。

窗口长度取值范围会受到信号长度和显示要求的限制，在允许范围内可以观察窗口长度对于分辨率的影响，具体分析见b.3。 窗函数类型 利用wavtool 函数可以查看各个窗口函数的图像

图v.2.b.1 Blackman 的窗函数图像如图

2.c 加窗 Fourie r 正变换案例分析

对比一下分析分时频信号

y ={

t <0.5, (sin (2?pi ?50?t1))t ≥0.5, (sin (2?pi ?120?t1))

图v.2.c.1 输入信号为一分时频信号，在t<0.5时，其表达式为y1=sin(2*pi*50*t1)

t>0.5时，其表达式sin(2*pi*120*t2)。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的等高线图。选取的窗函数为矩形窗

先分析右上角的ＦＦＴ幅度谱：

将源信号图放大后可观察到在采样点５１１处出现了相位的间断点，而在快速傅里叶变换中也可观测到幅度特性以此处为对称点。但是这不意味着快速傅里叶变换有很好的时间局部化特征。

图v.2.c.2放大的 ＦＦＴ幅度谱

图v.2.c.3放大的原函数

当输入一个分段三频率的信号时可见，时频信息被混在幅度谱里。 所以分时频信号不具备时频可分性，即无法从时域和频域中的任何一个域完整的描述信号的全部特征。

图v.2.d.1.输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor图。选取的窗函数为gausswin窗。窗长为20.

function timefreq2见代码code acht.

其次分析右下角的时频谱图,将其旋转放大后得到如下图:

图v.2.d.1.1 输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。此图为其时频谱图。

由其时间分辨率特性来看，在时间轴信号幅度不变。由其频域分辨率来看，其频率连续线性增长由0Hz~约60Hz.

可见wfＦＴ具备一定的时频特性。

窗长对于频率分辨率的影响：

图v.2.d.2.输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor图。选取的窗函数为gausswin窗。窗长为25.

当将窗长依次为20，25，30，40，60，80时（沿图2，3，1顺时针方向顺序），再次观察加窗傅里叶变换的时频谱图，发现分辨率越来越高，而且底色蓝越来越深，这意味着能量越少泄露，特性越好，在这里选择了最优的窗函数gausswin. 图v.2.d.3.输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。六幅图均为其时频谱图，图中横坐标为时间点n,纵轴为模拟频率Hz

分析窗函数的类型对时频分辨率的影响：

图v.2.c.6 输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor图。选取的窗函数为矩形窗。

观察可以看到等高图中的频率虽然线性增长，但是却有很多的抖动的干扰谱线频域分辨率不太好。

图v.2.c.7输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的等高线图。选取的窗函数为hamming窗。

相比较矩形窗，hamming窗的时频分辨率的更高，频率的线性增长更加清晰。而且对比左下角图，其带外更为平坦。

图v.2.c.8输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor 图。选取的窗函数为rectangular 窗。

图v.2.c.9.输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor 图。选取的窗函数为kaiser 窗。

图v.2.c.10.输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor 图。选取的窗函数为blackman 窗。

图v.2.c.11.输入信号为x = exp(j*4*pi*(t/80).^2)，复多普勒信号。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的pcolor 图。选取的窗函数为gaussian 窗。

图v.2.c.8（rec ） 和 图v.2.c.9(Kaiser) 都有非常明显的能量泄露。

图v.2.c.10 ( blackman) 和 图v.2.c.11（gaussian ）相比前两幅图，可看到能量泄露明显减少，但是付出了分辨率的代价。 其中gaussian 窗是最为折中的窗函数。 总结：

选择窗函数的的原则是： 1.窗谱的主瓣窄而高，以提高分辨率；

2.旁瓣幅值应该小，以减少频谱泄露。（1，2均指频谱特性）

当主瓣宽度较窄时，旁瓣幅度较高，能量泄露严重；但当旁瓣幅度小时，可以得到较平坦的带外响应但是主瓣宽度较宽。其具体选择应视具体情况而定，其根据为信号性质。

比如对于算法稳定性要求较高的电网电器设备自动检测的情况下，常选用Bartlett 窗和Hanning 窗，以较少初始相位的敏感性。

图v.2.c.12 图为矩形窗，Hamming 窗和Blackman 窗和gausswin 的时域和频域谱图。红色的为Blackman 窗，绿色的为Hamming 窗，蓝色的为rectwin 窗,绿色的为rectwin 窗。

矩形窗在时域内带内无衰减，带外衰减为零。在频域抖动剧烈，能量局部化程度较低。

Hamming 窗在时域内窗口宽度较宽，衰减较慢，在频域里带外衰减比矩形窗快。

Blackman 窗在时域内宽度比矩形窗和Hamming 窗都窄，率减速度很快。

在频域里主瓣比Hamming 窗窄，但是带外Blackman 窗的衰减程度比Hamming 快。

Gausswin 具有很好的时频局部化性质。 其他具体参数见表图v.2.c.13。

五：小波变换（Wavelet Transform ） 1.Wavelet transform 正变换数学描述

相比于加窗傅里叶变换，连续小波变换是一种自适应的时频窗结构。

将母小波伸缩平移，可以得到很多副本满足自适应的需求。其数学表达为：

ψa,τ(t )(t ?τ

a ),a >9，τ ? R 参数解释：ψa,τ(t )为小波基函数，a 为尺度因子和伸缩因子，τ为平移因子。 连续小波变换的定义：

将连续小波基函数ψa,τ(t )作用于能量有限信号f (t )：

WT f (a,τ)=?f (t ),ψa,τ(t )?a

f (t )ψ?(t ?τ

a )?t R

WT f (a,τ)为小波变换系数。

2.Wavelet transform 正变换Matlab 实现

可以使用matlab 的cwt 函数，以下内容使用的是matlab 自带的小波工具箱图fc.1.2.坐标分别为，scale ，b 和幅度值。

cwt 函数为可实现一维连续小波变换：

COEF = cwt( S, SCALES,‘wname‘,‘plotmode‘)

plotmode 可以通过的MATLAB 的 axe properties 调整着色模式。 代码见c ode fünf

离散化后的小波变换可写为：

ψm,n (x )=a 0?

m

2Ψ(a 0?m (x ?nb 0a 0m ))=a 0?

m

2Ψ(a 0?m

x ?nb 0) s

a 0?m

为伸缩因子，b 平移因子。

通过观察横坐标，横坐标为b,纵坐标为scale,a.a 越大，则伸缩因子越大，窗口越大，频率越低。b 越大,时间参数越大。

3.Wavelet transform 正变换分析总结

以exp(j*4*pi*(t/80).^2)为例

从图f.3..1中可看出信号的频率在随时间增长。

图f.3..2如图为实多普勒信号的连续小波变换图谱。Level=4

图f.3..3如图为实多普勒信号的连续小波变换图谱。Level=3

图f.3..4如图为实多普勒信号的连续小波变换图谱。Level=2

对比图fc.1.1~3可看出在分析多普勒信号时，在一定范围内，Level 越

大信号的频率分辨率越好。

一下选取Daubechies,complex gaussian, symlet 三种函数基本同尺度分解同一信号。结果如下：

图f.3..5如图为实多普勒信号的连续小波变换图谱。Level=4

wavefun=Daubechies

图f.3..6如图为实多普勒信号的连续小波变换图谱。Level=4

wavefun= Complex gaussian

图f.3..7如图为实多普勒信号的连续小波变换图谱。Level=4 wavefun= symlet 代码见

code sechs

对比可见在分析多普勒信号时，Complex gaussian的时间频率分辨率最佳。

在分析时，小波函数的选取依照信号的特征各异。

4. 实测语音等信号的案例分析

在此部分使用小波变换工具包，MATLAB工具包如下：

我们使用wavelet 1D,因为音频信息一般是二维的，如果是双声道则通过矩阵运算分开处理。

此音频文件被采样的矩阵尺度为48000的单声道矩阵，具体信息如下：

图为原音频文件的拟合信号，采用haar小波进行了尺度为5的分解。

图v.1.1为sampledataB（男生所说“老师新年快乐”）的分析结果。

图v.1.2为Full decomposition 显示模式下的五层分解结果。

图v.1.3为sampledataB的连续小波变换，小波函数为haar.具体信息如下：

图v.1.4,第一幅图为sampledataG的wavelet-1D harr五层分解，第二幅图为其多尺度分解的结果。第三幅图为连续小波分析的结果。

对比wavelet 1D 和continous wavelet

的信号分析都可以看到男孩的声音能

量较低，能量较高的频带也比女性声音主要频带低，但是每个音节里，男性的声音音域比女生宽广，即谱线低高频系数基本相等。这大概是为什么男性的声音低沉，而女性的声音比较灵动的原因。

图v.2.1信号erhu的九层haar函数分析结果。

图v.2.2信号paino的九层haar函数分析结果

图v.2.3信号hulusi的5层haar函数分析结果

信号来源：歌曲“风居住的街道”源音频文件通过Adobe Audition 剪辑，在经过矩阵运算提取了第一个声道，具体音频请听附件。

wavelet多尺度分析:

从细节信号（d）的电平变化趋势中：

首先在钢琴的细节信号中，电平呈现出近似指数衰减的形态，电平爬升很快，近乎冲激但是衰减同样也很快，所以近似指数衰减。

其次在葫芦丝的细节信号中，电平呈现出近似纺锤型的形态，爬升很慢，衰减也很慢。其中一个音看似衰减很快那是因为钢琴伴奏的缘故。

最后在二胡的细节信号中，电平呈现前两者结合的形态，但是电平爬升仍然比衰减要慢一点。

我们猜测这和乐器的发生机理有关系。

钢琴是按下键盘上的琴键，牵动钢琴里面包着绒毡的小木槌，继而敲击钢丝弦发出声音。所以电平升高很快，因为阻尼较大所以衰减也很快。

葫芦丝是靠着气流冲击簧片，簧片震动，并在管身腔体共鸣震动所产生的。所以人所吹的气体传感到簧片，簧片的跟随速度自然比钢琴直接接触的慢很多。而且空气柱震动在空气中衰减也比较慢，所以音色比较悠扬。

而二胡是由音箱共振音箱共震，由于按弦使琴弦震动的长度不同，震动的频率改变，使共鸣的音箱震动频率不同，由音箱放大发出声音! 而二胡显而易见在频率跳转时电平收敛速度较快。而因为共鸣箱震动使得电平爬升也比较慢。

图v.2.4信号paino的continous wavelet haar函数分析结果

图v.2.5信号hulusi（葫芦丝）的continous wavelet haar函数分析结果.

图v.2.5信号erhu的continous wavelet haar函数分析结果

Continous wavelet 分析

虽然没有条件对各个乐器进行单音频分析，但是由此依然可以看出个个乐器的音频特性。

首先，钢琴清晰可见每条谱线在低高频的分部近乎相等。也许这是所谓钢琴音域宽广的原因，我们猜测这也许和期冲激爬升的特性有关，因为冲激函数的频谱是白色的。

其次葫芦丝，主要的能量分布在低频区，高频区也有，但是中频似乎被陷波了。低频猜测和气体发生的激励有关。而高频的分量也许是葫芦丝的银色比较清丽的原因。

最后二胡则每个音符都有不小的低频分量但是并不知道哪里是它的音色特征。

所以我们可以选取音色十分相近的小提琴和二胡做进一步的分析，因为小提琴的银色被认为是优美婉转的，相比二胡则沉郁悲怆。

图v.2.6信号LZviolin 的continous wavelet haar函数5层分析结果

图v.2.7信号LZerhu 的continous wavelet haar函数5层分析结果

（信号来源，梁祝主题截取，刚好避开了所有的伴奏，相信请听附件。）在此仅对频谱进行分析：

图v.2.6可看到，小提琴的音域比二胡广一些，谱线从低频到高频分布较二胡均匀。也许这是他比较婉转的原因，比较葫芦丝钢琴和小提琴可以粗略估计音域宽广的乐器也许更加让人觉得婉转。

从图v.2.7可看到能量局部化的现象极其明显，各个音符集中在特定的时间和频率区域。而且高频与低频系数的差异比小提琴要大很多。而且高频的音符的局部化比低频要明显很多。这种音域局部化，也许是二胡听起来比较凄婉的原因。

参考文献：

[1]周伟.MATLAB小波分析高级技术.[M]西安电子科技大学出版社

[2]董长虹.Matlab小波分析工具箱原理与应用. [M]国防工业出版社

[3]张德丰.Matlab小波分析. [M]机械工业出版社

[4] Ingrid Daubechies. Ten Lectures on wavelets中文版[M]国防工业出版社

[5] Lokenach Debnach. Introduction to Hilbert Space [M] ELSEVIER ACADEMIC PRESS

附录：

Code erst

x=[......]; % discrete signal in time domain

N = length(x)

Mn = dftmtx(length(x)); %generate a dft matrix

Mn

length(Mn)

Fx =Mn*x;

subplot (1,2,1);

stem (x, 'fill');

subplot (1,2,2);

stem(abs(Fx), 'fill');

code zwei

clear all;

%% FS sampling rate

%T sampling time % L signal width % t time series

Fs = 1000;

T = 1/Fs;

L = 1000;

t = (0:L-1)*T;

%% input signal 50+120Hz sin

x = sin(2*pi*50*t)+sin(2*pi*120*t);

y = x ; %+ 2*rand(size(t)); % add a ramdom signal

subplot(1,2,1);

plot(Fs*t(1:50),y(1:50))

title('signal added with a random noise')

xlabel('time/ms');ylabel('y');

%% FFT

NFFT = 2^nextpow2(L);

Y = fft (y,NFFT)/L;

f = Fs/2*linspace(0,1,NFFT/2+1)

subplot(1,2,2);

plot(f,2*abs(Y(1:NFFT/2+1)))

title('y(t)')

xlabel('Hertz');

ylabel('Magnitude');

code drei

clc;clear all;

DF=100 t1=-DF/2-1: 0.1: 0; % t<0

y1=sin(2*pi*t1/5);

%figure(1)

%plot(t1,y1)

t2=0.1: 0.1: DF/2+1; % t>0

y2=sin(2*pi*t2/15);

% figure(2)

% plot(t2,y2)

t=[t1 t2] %-DF/2-1: 0.1: DF/2+1;

%t

y=[y1 y2] %·?ê±?μD?o? t<0ê± T=5,t>0ê± T=15

N=400;

x = zeros(1,N);

t = 0:N-1;

x = y; %signal to be analysed

figure(1);

timefreq(x,20,'rec');

figure(2);

timefreq(x,20,'Blackman');

funtion timefreq

%% funtion timefreq()

%step1- plot of the test signal

function timefreq(x,Nw,window)

% explain the parameters

% x test signal , a matrix

% Tw width of Time Window

% window type of the window that we may choose to fetch the time pieces

%% test signal 3 frequency

subplot(2,2,1);

plot (real(x));

title('real part of the test signal')

grid

X=fft(x);

X=fftshift(X); %μ÷??á??μ

subplot(2,2,2)

plot(abs(X));

title('FFT of the signal, amplitude respond ')

grid % amplitude respond

%% add rectangular window

Lap=Nw/2; % reused length

Tn=(length(x)-Lap)/(Nw-Lap); % fractional length

nfft=2^ceil(log2(Nw)); % Number of fft

TF=zeros(Tn,nfft); % line = length of the fractional length column = frequency point number

for i=1:Tn %cycle shift the time domain window function if (strcmp(window,'rec'))

xw=x((i-1)*10+1:i*10+10);

elseif (strcmp(window,'hamming'))

xw=x((i-1)*10+1:i*10+10).*hamming(Nw)';

elseif (strcmp(window,'blackman'))

xw=x((i-1)*10+1:i*10+10).*blackman(Nw)';

elseif (strcmp(window,' hann '))

xw=x((i-1)*10+1:i*10+10).*hann（Nw)';

elseif (strcmp(window,' Kaiser '))

xw=x((i-1)*10+1:i*10+10).*Kaiser(Nw)';

elseif (strcmp(window,' Gausswin'))

xw=x((i-1)*10+1:i*10+10).* Gausswin (Nw)';

else

return;

end

temp=fft(xw,nfft) % FFT transfrom of the signal limited in the window function in time domain

temp=fftshift(temp); %adjust te location

TF(i,:)=temp;

end

%% plot the result of TF

subplot (2,2,3);

fnew=((1:nfft)-nfft/2)/nfft;

tnew=(1:Tn)*Lap;

[F,T]=meshgrid(fnew,tnew); %draw a mesh of both frequency domain and time domain-3D

grid

mesh(T,F,abs(TF));

xlabel('n');ylabel('\omega');zlabel('G_f')

title('time frequency respond')

subplot(2,2,4)

contour(T,F,abs(TF));

grid

xlabel('n');ylabel('\omega');

title('contour of the time frequency respond')

ｃｏｄｅｄｒｅｉ

DF=100

t1=-DF/2-1: 0.1: -20; % t<0

y1=sin(2*pi*t1/5);

%figure(1)

%plot(t1,y1)

t2=-19.9: 0.1: 20; % t>0 %(from the result )Columns 312

-19.9000

y2=-sin(2*pi*t2/15);

% figure(2)

% plot(t2,y2)

t3=20.1: 0.1: 100;

y3=-cos(2*pi*t3/60); %Columns 711 through 712

20.0000 20.1000

% broken point cheque

t=[t1 t2 t3]; %column 1021

%t

y=[y1 y2 y3]; %·?ê±?μD?o? t<0ê± T=5,t>0ê± T=15 %figure(2),coefs = cwt (y,1:0.1:200,'haar','plot')

%figure(3),mesh(res);

N=400;

x = zeros(1,N);

% t = 0:N-1;

x = y; %signal to be analysed

figure(1);

timefreq(x,20,'rec');

figure(2);

timefreq(x,20,'Blackman');

code vier

clc;clear all;

N=400;

x = zeros(1,N);

t = 0:N-1;

x = exp(j*4*pi*(t/80).^2);

figure(1);

timefreq(x,20,'rec');

figure(2);

timefreq(x,20,'hamming');

figure(3);

timefreq(x,20,'blackman');

code fünf

t=0:0.01:50;

f = sin(j*4*pi*(t/80).^2);

%f=cos(j*4*pi*(t/80).^2);

% figure(4);

% plot(f);

figure (2),

res= cwt ( f, 1:0.1:200, 'db4', 'plot' );

code sechs

t=0:0.01:50;

f = sin(0.1*t.^2);

%f=cos(j*4*pi*(t/80).^2);

%figure(4);

% plot(f);

figure(1),res= cwt(f,1:0.1:200,'db4','plot');

figure(2),res= cwt(f,1:0.1:200,'cgau4','plot');

figure(3),res= cwt(f,1:0.1:200,'sym4','plot');

code sieben

%% step1 sampledata hold the audio signal as a matrix￡?FS is sampling rate,44100Hz

[sampledata,FS]=audioread('C:\Users\é?3?\Documents\matlab\hom ework\Elizabeth.mp3')

%% step2 if double track ,then choose one of them

% then run size(sampledata) = 6797424 *2

% cho_matrix = [1; 0]

cho_matrix = [1; 0];

% convert the sampledata into single track

sampledata = sampledata * cho_matrix

% 2*6797424 = 2*6797424 * 2*1

>> load gong.mat

sound(y,2*FS)

>> size(y)

ans =

42028 1

Code acht

%% funtion timefreq()

%step1- plot of the test signal

function timefreq(x,Nw,window)

% explain the parameters

% x test signal , a matrix

% Tw width of Time Window

% window type of the window that we may choose to fetch the time pieces

%% test signal 3 frequency

subplot(2,2,1);

plot (real(x));

title('real part of the test signal')

grid

X=fft(x);

X=fftshift(X); %μ÷??á??μ

subplot(2,2,2)

plot(abs(X));

title('FFT of the signal, amplitude respond ')

grid % amplitude respond

%% add rectangular window

Lap=Nw/2; % reused length

Tn=(length(x)-Lap)/(Nw-Lap) % fractional length

nfft=2^ceil(log2(Nw)) % Number of fft

TF=zeros(Tn,nfft); % line = length of the fractional length column = frequency point number

for i=1:Tn %cycle shift the time domain window function

if (strcmp(window,'rec'))

xw=x((i-1)*Nw/2+1:i*Nw/2+Nw/2); % 1-10 10 signal point each time

elseif (strcmp(window,'hamming'))

xw=x((i-1)*Nw/2+1:i*Nw/2+Nw/2).*hamming(Nw)' ;

elseif (strcmp(window,'blackman'))

xw=x((i-1)*Nw/2+1:i*Nw/2+Nw/2).*blackman(Nw)' ;

elseif (strcmp(window,'gausswin'))

xw=x((i-1)*Nw/2+1:i*Nw/2+Nw/2).*gausswin(Nw)' ;

elseif (strcmp(window,'kaiser'))

xw=x((i-1)*Nw/2+1:i*Nw/2+Nw/2).*kaiser(Nw)' ;

elseif (strcmp(window,'hann'))

xw=x((i-1)*Nw/2+1:i*Nw/2+Nw/2).*hann(Nw)' ;

else

return;

end

temp=fft(xw,nfft); % FFT transfrom of the signal limited in the window function in time domain

temp=fftshift(temp); %adjust te location

TF(i,:)=temp;

end

%% plot the result of TF

subplot (2,2,3);

fnew=((1:nfft)-nfft/2)*400/nfft ; %

plot(Fs*t(1:50),y(1:50))

tnew=(1:Tn)*Lap;

[F,T]=meshgrid(fnew,tnew); %draw a mesh of both frequency domain and time domain-3D

grid

mesh(T,F,abs(TF));

xlabel('n');ylabel('\omega/Hz');zlabel('G_f')

title('time frequency respond')

%contour

subplot(2,2,4)

%contour(T,F,abs(TF));

%grid

%xlabel('n');ylabel('\omega');

%title('contour of the time frequency respond')

pcolor(T,F,abs(TF));

grid

xlabel('n');ylabel('\omega');

title('pcolor of the time frequency respond')

code neun

clc;clear all;

N=400;

x = zeros(1,N);

t = 0:0.5:N-1;

x = exp(j*4*pi*(t/80).^2);

subplot(2,3,1)

timefreq(x,20,'gausswin');

grid

title('pcolor winlength=20')

subplot(2,3,2)

timefreq(x,25,'gausswin');

grid

title('pcolor winlength=25')

subplot(2,3,3)

timefreq(x,30,'gausswin');

grid

title('pcolor winlength=30')

subplot(2,3,4)

timefreq(x,40,'gausswin');

grid

title('pcolor winlength=40')

subplot(2,3,5)

timefreq(x,60,'gausswin');

grid

title('pcolor winlength=60')

subplot(2,3,6)

timefreq(x,80,'gausswin');

grid

title('pcolor winlength=80')

code last

%% step1 sampledata hold the audio signal as a matrix￡?FS is sampling rate,44100Hz

[erhu,FS]=audioread('C:\Users\é?3?\Documents\matlab\homework \?toú_·??ó×?μ???μà_???-.mp3')

[paino,FS]=audioread('C:\Users\é?3?\Documents\matlab\homework \???ù_·??ó×?μ???μà_???-.mp3')

[hulusi,FS]=audioread('C:\Users\é?3?\Documents\matlab\homewor k\où????_·??ó×?μ???μà_???-.mp3')

%% step2 if double track ,then choose one of them

% compare the erhu hulusi and paino

[erhu,FS]=audioread('C:\Users\é?3?\Documents\matlab\homework \?toú_·??ó×?μ???μà_???-.mp3')

cho_matrix = [1; 0];

% convert the sampledata into single track

erhu = erhu * cho_matrix

[paino,FS]=audioread('C:\Users\é?3?\Documents\matlab\homework \???ù_·??ó×?μ???μà_???-.mp3')

cho_matrix = [1; 0];

paino = paino * cho_matrix

[hulusi,FS]=audioread('C:\Users\é?3?\Documents\matlab\homewor k\où????_·??ó×?μ???μà_???-.mp3')

cho_matrix = [1; 0];

hulusi = hulusi * cho_matrix

%% compare the violin and erhu

[LZerhu,FS]=audioread('C:\Users\é?3?\Documents\matlab\homewo rk\erhuáo×￡???-°? .mp3')

cho_matrix = [1; 0];

LZerhu = LZerhu * cho_matrix

[LZviol,FS]=audioread('C:\Users\é?3?\Documents\matlab\homewor k\violináo×￡???-°? .mp3')

cho_matrix = [0; 1];

LZviol = LZviol* cho_matrix

%?C:\Users\é?3?\Documents\matlab\homework\erhuáo×￡???-°? .mp3

%?C:\Users\é?3?\Documents\matlab\homework\violináo×￡???-°? .mp3