山东交通学院毕业设计(论文)
谱却是是一样的，这是因为傅里叶变换是一种整体变换，即对信号的表征要么完全在时域，要么完全在频域，作为频域表示的功率谱并不能告诉我们其中某种频率分量出现在什么时候及其变化情况。然而，在许多实际应用场合，信号是非平稳的，其统计量(如相关函数、功率谱等)是时变函数。这时，只了解信号在时域或频域的全局特性是远远不够的，最希望得到的乃是信号频谱随时间变化的情况。为此，需要使用时间和频率的联合函数来表示信号，这种表示简称为信号的时频表示。时频分析的主要研究对象是非平稳信号或时变信号，主要的任务是描述信号的频谱含量是怎样随时间变化的。时频分析有短时傅里叶变换(Short-time Fourier transform，简记为STFT)、Gabor展开和小波变换(Wavelet Transform，简记为WT)等。本文主要应用的方法是短时傅里叶变换。
为了分析语音信号，Koenig等人提出了语谱图(Spectrogram)方法2，定义为信号的短时傅立叶变换STFT的模平方，故亦称为STFT方法或者STFT谱图[9]。离散短时傅立叶变换定义如下：
式中
是时间窗函数(一般为汉明窗)。短时傅立叶变换的基本思想是用一
个时间宽度足够窄的固定的窗函数乘时间信号，使取出的信号可以被看成平稳的，然后对取出的一段信号进行傅立叶变换，便可以反映出该时间宽度中的频谱变化规律，如果让这个固定的窗函数沿着时间轴移动，那就可以得到信号频谱随时间变化的规律了。
下面做出上述两个已知信号的时频图。 程序4：
N=11025; %采样点数 f1=100;f2=200; fs=11025;
t=0:1/fs:((N-1)*1/fs);
x=sin(2*pi*f1*t)+sin(2*pi*f2*t); subplot(2,1,1);plot(t,x);grid on; xlabel('时间(秒)');ylabel('幅度(伏)'); xlim([0 ((N-1)*1/fs)]); title('时域波形图'); window=512;%窗长度 noverlap=window/2;%重合度 2 语谱图即时频图，这里spectrogram函数同样适用于音乐信号。
8
山东交通学院毕业设计(论文)
F=window;
[S,Fd,Td,P] = spectrogram(x, window, noverlap, F ,fs); subplot(2,1,2);
surf(Td,Fd,20*log10(P),'edgecolor','none'); axis tight; view(0,90); xlabel('时间 t(s)'); ylabel('频率 f(Hz)'); ylim([0 fs/20]) title('FFT变换时频图'); set(gca,'YDir','normal'); colorbar; 结果如下：
时域波形图21幅度(伏)0-1-200.10.20.30.50.6时间(秒)FFT变换时频图0.40.70.80.9 频率 f(Hz)400-100200-2000 0.10.20.30.40.50.6时间 t(s)0.70.80.9-300
图1.4信号1的波形和时频图
Fig. 1.4 The waveform and time-frequency diagram of signal 1
同样我们绘出信号2的时频图(程序见附录)，结果如下：
9
山东交通学院毕业设计(论文)
时域波形图10.5幅度(伏)0-0.5-100.20.40.611.2时间(秒)FFT变换时频图0.81.41.61.8 -50-100-150-200频率 f(Hz)4002000 0.20.40.60.811.2时间 t(s)1.41.61.8
图1.5信号2的波形和时频图
Fig. 1.5 The waveform and time-frequency diagram of signal 2
需要指出的是：第一，长窗具有较高的频率分辨率，但具有较低的时间分辨率。从一个周期到另一个周期，共振峰是要发生变化的，这一点即使从信号波形上也能够看出来。然而，如果采用较长的窗，这种变化就模糊了，因为长窗起到了时间上的平均作用。第二，短窗的频率分辨率低，但具有较高的时间分辨率。采用短窗时，能够从短时频谱中提取出共振峰从一个周期到另一个周期所发生的变化。当然，激励源的谐波结构也从短时频谱上消失了。所以，在对信号进行短时傅里叶分析时，需要折衷考虑时间分辨率与频率分辨率以适当选取窗长[10]。
下面我们做一个简单的时频分布图，并取不同的窗长来说明 程序5：
filename='E:\\音频文件\\piano-G4.wav'; [x,fs]=wavread(filename); t=0:1/fs:((length(x)-1)*1/fs); N=length(x); NFFT = 4096 Y = fft(x,NFFT)/N;
f = fs/2*linspace(0,1,NFFT/2+1); subplot(3,1,1);
10
山东交通学院毕业设计(论文)
plot(f,2*(abs(Y(1:NFFT/2+1)))); xlim([0 fs/8]); xlabel({'Frequency'}) grid on;
title('FFT变换频谱图'); clear f;clear Y;
window=2048;%窗长度 noverlap=window/2;%重合度 F=window;
[S,Fd,Td,P] = spectrogram(x, window, noverlap, subplot(3,1,2);
surf(Td,Fd,20*log10(P),'edgecolor','none'); axis tight; view(0,90); xlabel('时间 t(s)'); ylabel('频率 f(Hz)'); ylim([0 fs/12])
title('FFT变换时频图-窗长2048'); set(gca,'YDir','normal'); colorbar;
window=256;%窗长度 noverlap=window/2;%重合度 F=window;
[S,Fd,Td,P] = spectrogram(x, window, noverlap, subplot(3,1,3);
surf(Td,Fd,20*log10(P),'edgecolor','none'); axis tight; view(0,90); xlabel('时间 t(s)'); ylabel('频率 f(Hz)'); ylim([0 fs/12])
title('FFT变换时频图-窗长256'); set(gca,'YDir','normal'); colorbar; 运行结果如下：
11
F ,fs); F ,fs);