信息工程概论窗口傅里叶变换.docx

信息工程概论窗口傅里叶变换.docx

《信息工程概论窗口傅里叶变换.docx》由会员分享，可在线阅读，更多相关《信息工程概论窗口傅里叶变换.docx（10页珍藏版）》请在冰豆网上搜索。

信息工程概论窗口傅里叶变换

信息工程概论作业

——窗口傅里叶变换

姓名：

白子轩学号：

2130602008班级：

信计31

1、传统的傅里叶变换

我们都知道，信号分析中最重要的两个参数是时间和频率，而我们一般所得到的信号表示形式都是

的形式，而我们可以通过传统的傅里叶变换

，可以把信号变为频域表示

。

但是，传统的傅里叶变换只对平稳的信号有用。

对于非平稳的信号需要用时间和品率的联合函数来表示信号。

因此，我们需要短时傅里叶变换，也就是窗口傅里叶变换。

2、窗口傅里叶变换

对于信号的频率是随时间变化的信号。

为了获得它的随时间变化的频谱，最采用的处理办法是加窗技术对信号截取，然后对截取的局部信号作Fourier变换。

然后不断地移动窗口函数中心的位置，就可以得到信号的局部区域的瞬时频率，

因此，对于连续的信号，它的窗口傅里叶变换为：

窗口傅里叶逆变换为：

而对于离散的信号，它的窗口傅里叶变换为：

窗口傅里叶逆变换为：

3、窗口函数

要进行窗口傅里叶变换，首先要要选择窗口函数，窗口函数有很多，例如高斯窗、hamming窗和Hanning窗等等。

其中高斯窗函数被设计为了分析瞬态信号，Hamming和Hann窗函数被设计为了分析窄带信号，Kaiser-Bessel窗函数可用于更好地分离两个频率成分非常接近但振幅完全不同的信号。

在matlab中我们也可以直接调用一些窗口函数，调用的方法如下：

4、实验

实验1：

题目：

在这里我先做了一下书上的例子，对线性调频信号

进行频谱分析。

题目分析：

这个例子有两种做法，第一种方法是直接调用matlab中的spectrogram函数，第二种方法是按照定义选取窗口函数，然后对每一小段的做快速傅里叶变换就可以了。

方法一：

源程序：

clear

t=0:

0.001:

10;

t1=t;

f1=sin（2*pi*2*power（t,2））;

subplot（2,2,1）;

plot（t,f1）;

subplot（2,2,2）;

g=1/6*exp（-0.5*power（t,2））.*（t>=-3&t<=3）+0.*（t>3|t<-3）;

t=-4:

0.01:

4;

g1=1/6*exp（-0.5*power（t,2））.*（t>=-3&t<=3）+0.*（t>3|t<-3）;

plot（t,g1）;

subplot（2,2,3）;

[S,F,T,P]=spectrogram（f1,gausswin（600）,580,600,1E3）;

surf（T,F,10*log10（abs（P））,'edgecolor','none'）;

axis（[010050]）

view（0,90）;

xlabel（'Time（Seconds）'）;

ylabel（'Hz'）;

subplot（2,2,4）;

surf（T,F,10*log10（abs（P））+80,'edgecolor','none'）;

axis（[0100500200]）

%axistight

xlabel（'Time（Seconds）'）;

ylabel（'Hz'）;

zlabel（'enargy'）

得到结果：

结果分析：

我们用的是高斯窗口，得到了一个很好的结果，无论是2D图还是3D图，都与书上的图十分相似。

但有一个十分大的缺陷，就是无法重构原来的信号，因为我们是直接调用的spectrogram函数，并不太知道里面的具体程序是什么样的，所以无法还原原信号，也无法计算误差。

因此我们就需要第二种方法。

方法二：

源程序：

clear

t=0:

0.001:

10;

t1=t;

f1=sin（2*pi*2*power（t,2））;

subplot（3,2,1）;

plot（t,f1）;

subplot（3,2,2）;

g=1/6*exp（-0.5*power（t,2））.*（t>=-3&t<=3）+0.*（t>3|t<-3）;

t=-4:

0.01:

4;

g1=1/6*exp（-0.5*power（t,2））.*（t>=-3&t<=3）+0.*（t>3|t<-3）;

plot（t,g1）;

N=length（f1）;

Nw=20;%窗函数长windowlength

L=19;%窗函数每次移动的样点数,重叠宽度

Ts=round（（N-Nw）/L）+1;%计算把数据x共分成多少段

nfft=2^ceil（log2（Nw））;%FFT的长度

TF=zeros（Ts,nfft）;%将存放三维谱图，先清零

%fori=1:

Ts

i=0;

flag=0;

whileflag==0;

i=i+1;

if（i-1）*Nw+Nw

%y（（i-1）*L+1:

i*L+L）

%hamming（Nw）

xw=f1（（i-1）*L+1:

i*L+L）*hamming（2*L）;%取一段数据

temp=fft（xw,nfft）;%FFT变换

temp=fftshift（abs（temp））;%频谱以0频为中心

%TF（i,:

）;

TF（i,:

）=temp;%把谱图存放在TF中

else

flag=1;

end

end

subplot（3,2,3）;

%mesh（abs（TF））;

%view（0,90）;

%axistight

%imagesc（TF）;

contour（abs（TF））;

xlabel（'时间'）;

ylabel（'频率'）

subplot（3,2,4）;

mesh（abs（TF））;%三维绘图

axistight

title（'STFT'）;

xlabel（'时间'）;

ylabel（'频率'）;

%%%%短时傅里叶变换

X=fft（f1）;

X=fftshift（X）;

%%%%重构

x1=ifftshift（X）;

x1=ifft（x1）;

subplot（3,2,5）;

plot（t1,x1）;

xlabel（'时间'）;%x轴

ylabel（'振幅'）;%y轴

%%%%误差

e=x1-f1;

subplot（3,2,6）;

plot（t1,abs（e））;

xlabel（'时间'）;%x轴

ylabel（'振幅'）;%y轴

得到结果：

结果分析：

我们发现得到的结果非常不理想，经过多次试验，我并没有找到是哪里出的问题，因此在接下来的实验中将要放弃这种方法。

实验2：

题目：

录一段自己的声音，对这段声音进行频谱分析。

实验分析：

这个题目和实验1类似，只不过一个用的是连续的窗口傅里叶变换，一个用的是离散的窗口傅里叶变换，所以在这里我们仍然运用matlab中的spectrogram函数即可。

在这里我们需要注意一下，mp3是双音轨的，所以数字化后是一个是一个二维数组，而我们只需要分析其中的一组即可。

于是我就录了自己的一段笑声，并抽取其中的一组数组，进行分析。

源程序：

clear

[x,fs]=audioread（'C:

\Users\lenovo\Desktop\小波分析\笑.mp3'）;

y=x（:

1）';

subplot（3,2,1）;

plot（y）;

title（'信号波形图'）;%图名

xlabel（'时间'）;%x轴

ylabel（'振幅'）;%y轴

holdon

N=length（y）;

Nw=600;%窗函数长windowlength

L=Nw/2;%窗函数每次移动的样点数,重叠宽度

%%%%%%窗口函数

subplot（3,2,2）;

g=hamming（Nw）;

j=-Nw/2;

fori=1:

Nw

t（i）=j;

j=j+1;

end

plot（t,g,'.'）;

subplot（3,2,3）;

[S,F,T,P]=spectrogram（y,hamming（600）,580,600,1E3）;

surf（T,F,10*log10（P）,'edgecolor','none'）;

axistight

view（0,90）;

xlabel（'Time（Seconds）'）;

ylabel（'Hz'）;

subplot（3,2,4）;

surf（T,F,10*log10（P）,'edgecolor','none'）;

axistight

xlabel（'Time（Seconds）'）;

ylabel（'Hz'）;

zlabel（'enargy'）

Y=fft（y）;

X=fftshift（Y）;

%%%%重构

x1=ifftshift（X）;

x1=ifft（x1）;

subplot（3,2,5）;

plot（x1）;

xlabel（'时间'）;%x轴

ylabel（'振幅'）;%y轴

%%%%误差

e=x1-y;

subplot（3,2,6）;

plot（abs（e））;

xlabel（'时间'）;%x轴

ylabel（'振幅'）;%y轴

得到结果：

结果分析：

我们可以看到得到了一个比较比较清晰地频谱图，但是它的3D图并不是太清晰。

结果不是太理想。