1、科研训练报告 一：科研训练目标及意义 在深入理解经典 Fourier 变换的基础上，通过系统介绍加窗 Fourier 变换（ Windowed Fourier Transfrom, WFT） ,小波变换（ Wavelet trsnaform,WT）等定义，理解时间分辨率，频率分辨率及，时频局部化性质等基本概念；通过编程实现基本的时频分析方法，理解时变信号中时频结构信息。 二：科研训练任务 1.理解加窗 Fourier 变换的意义 2.编程实现加窗 Fourier 分析几种合成信号 3.利用加窗 Fourier(小波变换 )分析几种 典型的合成信号 4.观察分析窗函数，窗长变化，频率，分辨率的影
2、响 5.可尝试利用加窗 Fourier（小波变换）正变换分析实测语音等信号 三 :报告结构 四：加窗 Fourier 分析 v.1.Fourier 正变换 v.1.a 离散 Fourier 正变换数学描述 有限长序列通过离散傅里叶变换，并通过快速傅里叶变换算法 FFT 即得：() = ()(1)(1)=1其中： = 2/ v.1.b.Fourier 正变换 Matlab 实现 被 分析信号为 N 维列向量； x = x(0),x(1),x(2),x(3)x(n1) 变换因子为： 0 0 0 (1)(1) DFT 结果为 X = (0),X(1),X(2),X(3)X(n1) 以矩阵形式表示 D
3、FT： = （具体代码参见代码 code erst） v.1.c.Fourier 正变案例分析总结 step1 检测代码的正确性 当输入一 个离散的冲激函数的 Delta 函数，其 DFT 变换的幅度谱为一等幅，幅度为 1 的白色谱。（参见图 a.3.1） Step2 分析信号特性 当输入一个正弦函数时，运行结果如下：幅度谱为一对称的离散冲激信号，冲激的横坐标为其频率，纵坐标为其幅度。由观测可看出在主频的附近有明显的频谱泄露。（参见图 a.3.2） 当输入等幅的正弦波叠加则运行结果显示，两个正弦波的幅值并不相同，存在明显的栅栏现象。（参见图 a.3.3） 当输入一分时频信号时，该信号表达式为：
4、 = 0.5时，其表达式 sin(2*pi*120*t2)。横坐标为模拟频率，纵坐标为幅值。 横坐标为模拟频率 ， 纵坐标为其幅度值。 2.加窗 Fourier 分析 2.a 加窗 Fourier 正变换数学描述 通过参见图 a.3.1 和图 a.3.3 和图 a.3.4 对比分析可见，傅里叶变换是一种针对 -(,)的全局变换。 而加窗 Fourier 变换的基本思想是在时域内进行分段，用位置不同的窗函数与原信号相乘。在时域内，时窗函数一般选取具有能量局部化的函数，选定一个基本窗函数后，通过将窗函数中心沿时间轴平移，得到一组窗函数，使得窗函数很好的覆盖时间轴。因此，时间窗等效提取了源信号的不同
5、时间内的信号而时间带外的信号衰减为零。 加窗 Fourier 变换认为在窗内信号是准平稳的，可采用平稳信号的分析方法，如频谱分析。 最早的 Gabor 的定义如下 ; 设， () 2()即 0 0.5时，其表达式 sin(2*pi*120*t2)。右上角图为快速傅里叶变换，左下为窗口变换的时频谱图，右下角为时频谱的等高线图。选取的窗函数为矩形窗 先分析右上角的幅度谱： 将源信号图放大后可观察到在采样点处出现了相位的间断点，而在快速傅里叶变换中也可观测到幅度特性以此处为对称点。但是这不意味 着快速傅里叶变换有很好的时间局部化特征。 图 v.2.c.2放大的 幅度谱 图 v.2.c.3放大的原函数
6、 当输入一个分段三频率的信号时可见，时频信息被混在幅度谱里。 所以分时频信号不具备时频可分性，即无法从时域和频域中的任何一个域完整的描述信号的全部特征。 图 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)，
7、复多普勒信号。 此图为其时频谱图 。 由其 时间分辨率特性来看， 在时间轴信号幅度不变 。 由其频域分辨率来看，其频率 连续 线性增长由 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顺时针方向顺序），再次观察加窗傅里叶变换的时频谱图，发现分
8、辨率越来越高，而且底色蓝越来越深，这意味着能量越少泄露，特性越好，在这里选择了最优的窗函数 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图。选取的窗函数为矩形窗。 观察可以看到等高图中的频率虽然线性增长，但是却有很多的抖动的干扰谱线频域分辨率
9、不太好。 图 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.输入信号为
10、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图 。选
11、取的窗函数为 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均指频谱特性） 当主瓣宽度较窄时，旁瓣幅度较高，能量泄露严重；但当 旁瓣幅度小时，可以得到较平坦的带外响应 但是主瓣宽度较宽 。 其具
12、体选择应视具体情况而定，其根据为信号性质。 比如 对于算法稳定性要求较高的 电网电器设备自动检测的情况下，常选用 Bartlett窗和 Hanning窗，以较少初始相位 的 敏感性 。 图 v.2.c.12 图为矩形窗， Hamming窗和 Blackman窗 和 gausswin的时域和频域谱图。红色的为 Blackman窗 ，绿色的为 Hamming窗，蓝色的为 rectwin窗 ,绿 色的为 rectwin窗 。 矩形窗在时域内带内无衰减，带外衰减为零。在频域抖动剧烈，能量局部化程度较低。 Hamming窗在时 域内窗口宽度较宽，衰减较慢，在频域里带外衰减比矩形窗快。 Blackman窗
13、在时域内宽度比矩形窗和 Hamming窗都窄，率减速度很快。在频域里主瓣比 Hamming窗窄，但是带外 Blackman窗的衰减程度比 Hamming快。 Gausswin具有很好的时频局部化性质。 其他 具体参数见表 图 v.2.c.13。 名称 最大旁瓣幅度 /dB 主瓣近似宽度 (/(M-1) 倍频程下降率 /(dB/OCT) Rectangular -13.3 4 -6 hann -31.5 8 -18 hamming -42.7 8 -6 blackman -58.1 12 -18 五 ： 小波变换 （ Wavelet Transform ） 1.Wavelet transform
14、 正变换数学描述 相比于加窗傅里叶变换，连续小波变换是一种自适应的时频窗结构。 将母小波伸缩平移，可以得到很多副本满足自适应的需求。其数学表达为： ,() = 1( ), 9， 参数解释： ,()为小波基函数， 为尺度因子和伸缩因子 ， 为平移因子 。 连 续小波变换的定义： 将连续小波基函数 ,()作用于能量有限信号 ()： (,) = (),() = 1()( )(,)为小波变换系数 。 2.Wavelet transform 正变换 Matlab 实现 可以使用 matlab 的 cwt 函数 ，以下内容使用的是 matlab 自带的小波工具箱 图 fc.1.2.坐标分别为， scale
15、， b 和幅度值。 cwt 函数为可实现一维 连续 小波变换： COEF = cwt( S, SCALES,wname,plotmode) plotmode 可以通过的 MATLAB 的 axe properties 调整着色模式。 代码见 code fnf 离散化后的小波变换可写为 ： ,() = 02(0(00) = 02(0 0) s 0为伸 缩因子 ， 平移因子 。 通过观察横坐标，横坐标为 b,纵坐标为 scale,a.a 越大，则伸缩因子越大，窗口越大，频率越低。 b 越大 ,时间参数越大。 3.Wavelet transform 正变换分析总结 以 exp(j*4*pi*(t/8
16、0).2)为例 从 图 f.31 中可看出信号的频率在随时间增长。 图 f.32 如图为实多普勒信号的连续小波变换图谱 。 Level=4 图 f.33 如图为实多普勒信号的连续小波变换图谱 。 Level=3图 f.34 如图为实多普勒信号的连续小波变换图谱 。 Level=2 对比 图 fc.1.13 可看出在分析多普勒信号时， 在一定范围内， Level 越大信号的频率分辨率越好。 一下选取 Daubechies,complex gaussian, symlet 三种函数基本同尺度分解同一信号。结果如下： 图 f.35 如图为实多普勒信号的连续小波变换图谱 。 Level=4 wavef
17、un=Daubechies 图 f.36 如图为实多普勒信号的连续小波变换图谱 。 Level=4 wavefun= Complex gaussian 图 f.37 如图为实多普勒信号 的连续小波变换图谱 。 Level=4 wavefun= symlet 代码见 code sechs 对比可见在 分析多普勒信号时， Complex gaussian 的时间频率分辨率最佳。 在分析时，小波函数的选取依照信号的特征各异。 4. 实测语音等信号的案例分析 在此部分使用小波变换工具包， MATLAB工具包如下： 我们使用 wavelet 1D,因为音频信息一般是二维的，如果是双声道则通过矩阵运算分开
18、处理。 此音频文件被采样的矩阵尺度为 48000 的单声道矩阵，具体信息如下： 图为原音频文 件的拟合信号，采用 haar 小波进行了尺度为 5 的分解。 图 v.1.1 为 sampledataB（ 男生所说“老师新年快乐” ） 的分析结果。 图 v.1.2 为 Full decomposition 显示模式下的五层分解结果。 图 v.1.3 为 sampledataB 的连续小波变换 ，小波函数为 haar.具体信息如下： 图 v.1.4,第一幅图为 sampledataG 的 wavelet-1D harr 五层分解，第二幅图为其多尺度分解的结果。第三幅图为连续小波分析的结果 。 对比
19、wavelet 1D 和 continous wavelet 的信号分析都可以看到男孩的声音能量较低，能量较高的频带也比女性声音主要频带低， 但是每个音节 里， 男性的声音 音域 比 女生 宽 广， 即 谱线 低高频 系数 基本相等。 这大概是为什么男性的声音低沉，而女性的声音比较灵动的原因。 图 v.2.1 信号 erhu 的九层 haar 函数分析结果。 图 v.2.2 信号 paino 的九层 haar 函数分析结果 图 v.2.3 信号 hulusi 的 5 层 haar 函数分析结果 信号来源： 歌曲 “ 风居住的街道 ” 源音频文件 通过 Adobe Audition 剪辑，在经过
20、 矩阵 运算 提取了 第 一个 声道 ， 具体 音频 请听 附件。 wavelet 多尺度分析 : 从细节 信号 （ d） 的 电平 变 化 趋势 中 ： 首先 在 钢琴的 细节信号 中， 电 平 呈现 出 近似 指数衰减的 形态， 电平 爬升很快， 近乎 冲 激 但是 衰减 同样也很快， 所以近似 指数 衰减。 其次 在 葫芦丝 的 细节信号中， 电平 呈现出 近似 纺锤型的 形态， 爬升 很慢， 衰减也很慢。 其中一个 音 看似 衰减 很快 那是因为 钢琴伴奏的缘故。 最后 在 二胡 的 细节信号中， 电平 呈现 前两者 结合 的 形态， 但是 电平 爬升 仍然比 衰减 要 慢一点。 我们
21、猜测 这和 乐器的 发生机理 有关系 。 钢琴 是 按下键盘上的琴键，牵动钢琴里面包着绒毡的小木槌，继而敲击钢丝弦发出声音。 所以 电平 升高 很快， 因为 阻尼 较大 所以衰减 也很快。 葫芦丝 是靠着气流冲击簧片，簧片震动，并在管身腔体共鸣震动所产生的 。 所以 人所吹的气体 传 感到 簧片， 簧片 的 跟随速度 自然比 钢琴直接接触 的慢很多 。 而且 空气 柱 震动 在空气中 衰减 也比较慢 ， 所以 音色 比较悠扬。 而 二胡 是由 音箱共振 音箱共 震，由于按弦使琴弦震动的长度不同，震动的频率改变，使共鸣的音箱震动频率不同，由音箱放大发出声音 ! 而 二胡 显而易见 在频率 跳转
22、时 电平 收敛 速度 较快。 而因为 共鸣箱 震动 使得 电平 爬升 也比较慢 。 图 v.2.4 信号 paino 的 continous wavelet haar 函数分析结果 图 v.2.5 信号 hulusi（ 葫芦丝 ） 的 continous wavelet haar 函数分析结果 . 图 v.2.5 信号 erhu 的 continous wavelet haar 函数分析结果 Continous wavelet 分析 虽 然没有条件 对各个 乐器 进行 单音 频 分析 ， 但是 由此 依然可以看出个个 乐器 的 音频特性。 首 先， 钢琴 清晰可见 每条 谱线 在 低高频 的
23、分部近乎相等。 也许这是所谓 钢琴 音域宽广的 原因， 我们猜测 这也许和 期 冲激 爬升 的 特性有关，因为 冲激函数 的 频谱 是白色的。 其次 葫芦丝， 主要的 能量 分布在 低频 区， 高频 区 也有， 但是 中频 似乎 被 陷 波 了。 低频猜测 和 气体 发生的 激励有关。 而高频 的分量也许是 葫芦丝 的 银色 比较 清丽 的原因。 最后 二胡 则每个音符 都有 不小 的 低频分量 但是 并不知道 哪里 是 它的 音色 特征 。 所以 我们可以 选取 音色 十分相近 的 小提琴 和 二胡 做进一步的分析，因为 小提琴的 银色 被认为 是 优美婉转 的， 相比 二胡 则 沉郁 悲怆
24、。 图 v.2.6 信号 LZviolin 的 continous wavelet haar 函数 5 层 分析结果 图 v.2.7 信号 LZerhu 的 continous wavelet haar 函数 5 层 分析 结果 （ 信号来源， 梁祝 主题 截取， 刚好避开 了 所有 的 伴奏 ， 相信 请听 附件 。 ） 在此 仅对 频谱进行分析： 图 v.2.6 可看到， 小提琴 的 音域 比 二胡 广 一些， 谱线 从 低频 到高频分布较 二胡 均匀 。 也许这是 他比较婉转的原 因， 比较葫芦丝 钢琴 和小提琴 可以粗略 估计 音域 宽广的 乐器 也许更加 让人觉得 婉转。 从 图 v
25、.2.7 可看到 能量局部化的 现象 极其明显， 各个音符 集中在 特定的 时间和 频率 区域。 而且 高频 与 低频 系数的 差异 比 小提琴要 大 很多。而且 高频的 音符 的 局部化 比 低频 要明显 很多。 这种 音域 局部化， 也许是二胡 听起来 比较 凄婉 的原因 。 参考文献： 1周伟 .MATLAB小波分析高级技术 .M西安电子科技大学 出版社 2董 长虹 .Matlab小波分析工具箱 原理与应用 . M国防 工业 出版社 3张德 丰 .Matlab小波分析 . M机械工业出版社 4 Ingrid Daubechies. Ten Lectures on wavelets中文版
26、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
27、 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 w
28、ith a random noise) xlabel(time/ms);ylabel(y); % FFT NFFT = 2nextpow2(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; % t0 y2=sin(2*pi*t2/15); % figure(
29、2) % plot(t2,y2) t=t1 t2 %-DF/2-1: 0.1: DF/2+1; %t y=y1 y2 % t0 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)
30、 % 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(
31、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=2ceil(log2(Nw); % Number of fft TF=zeros(Tn,nfft); % line = length of the fractional length column = frequency point numb
32、er 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)
33、.*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 l
34、ocation 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) co
35、ntour(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; % t0 %(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 % br
36、oken point cheque t=t1 t2 t3; %column 1021 %t y=y1 y2 y3; % t0 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
37、; 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 fnf 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
38、 ); 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 matrixFS is samplin
39、g rate,44100Hz sampledata,FS=audioread(C:UsersDocumentsmatlabhomeworkElizabeth.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*679
40、7424 = 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
41、 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=(lengt
42、h(x)-Lap)/(Nw-Lap) % fractional length nfft=2ceil(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
43、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
44、).*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); fn
45、ew=(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
46、%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
47、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(
