目录
一、信号的傅里叶变换频谱分析
MATLAB中可以用快速Fourier变换(Fast Fourier Transformation,FFT)实现频谱求解。例如,对于一个简单分段线性函数,对其频谱进行分析:
程序如下:
fun = @(x)0.*(x<0 |x>5)+x.*(x>=0 & x<1)+1.*(x>=1 & x<4) +(5-x).*(x>=4 & x<5);%分段函数
x=-5:0.2:10;%x自变量
y=fun(x);%得到输入的平稳信号y
Dfy = fft(y);%离散Fourier变换
Dfy_shift = fftshift(Dfy); %对称变换得到对称的Fourier频谱
figure
subplot(1,3,1),plot(y); axis([0,80,-0.05,1.05]),xlabel('n'),ylabel('y');title('(a)原始平稳信号')
subplot(1,3,2),plot(abs(Dfy)),xlabel('n'),ylabel('|Dfy|');title('(b)离散频谱幅值信息')
subplot(1,3,3),plot(abs(Dfy_shift)),xlabel('n'),ylabel('|Dfy|');;title('(c)对称变换后的频谱')
![](https://i-blog.csdnimg.cn/blog_migrate/962f9256201a27144c7d04882a270f83.png)
可以看出,一个简单的平稳信号如图1(a)所示,其频谱图很丰富,它的Fourier频谱如图1(b)所示,对称变换后的频谱如图1(c)所示。要想知道某一点的频率,需要利用原始信号所有点的值来得到。信号在某时刻的小区间上发生变化,信号的整个频谱都要受到影响,而频谱的变化无法标定发生的时间位置和变化的剧烈程度。这也就是Fourier变换的整体性之所在,即对信号的奇性不敏感。为改变这一点,1946年D.Gabor提出的窗口Fourier变换,即Gabor变换,实现了时频同时分析,给出了局部时间内的频谱信息的描述。
二、Gabor变换——窗口傅里叶变换
Gabor变换的性质分析:
·变换的窗口的形状大小是一定的,不随ω的变化而变化。
·变换离散化后不再正交。
·变换也有相应的重构公式
Gabor变换的程序实现如下:
- 选择合适的窗口函数,确定控制窗口的参数。
- 由原信号与窗口函数的乘积得到新的函数。
- 用FFT对新函数做Fourier变换。
- 得到对应的Gabor变换频谱。
例如,对稳信号的Gabor变换频谱。
现有一个图2(a)所示的简单平稳信号,其 Gabor 变换频谱和对称变换后的频谱的立体图如图2(b)和2 (d)所示。对应频谱的图像形式如图2(c)和2-5(e)所示。从 Fourier 频谱图(图1(c))和 Gabor 平面频谱可以看出,Gabor 变换的频谱更加集中,实现了时频的共同分析。但是该分析仍然受窗口函数的影响,时间和频率的分析相互独立,还没有从根本上突破Fourier 分析的桎梏。小波分析将改变这一点。
![](https://i-blog.csdnimg.cn/blog_migrate/499c1ef7c2880118ce72425fc6993e43.png)
matlab程序
x=-5:0.2:10;%x自变量
a = 1; %Gauss窗口函数的窗口大小参数
b = -5:0.2:10; %平移参数,一般为原始信号长度。
Dfy = zeros(length(b),length(b));
fun = @(x)0.*(x<0 |x>5)+x.*(x>=0 & x<1)+1.*(x>=1 & x<4) +(5-x).*(x>=4 & x<5);%分段函数
y=fun(x); %得到输入的平稳信号y
for i = 1:length(b)
GDfy(i,:) = fft(Gb_fun(a,b(i),x));%针对不同的b的离散Fourier变换
GDfy_shift(i,:) = fftshift(GDfy(i,:));%平移频率中心
end
%绘制频谱立体图和图像
figure
mesh(abs(GDfy)); xlabel('b'),ylabel('w'),zlabel('|Gabor(a,w)|');
figure
imshow(abs(GDfy));
%绘制对称变换后的频谱立体图和图像
figure
mesh(abs(GDfy_shift)); xlabel('b'),ylabel('w'),zlabel('|Gabor(a,w)|');
figure
imshow(abs(GDfy_shift));
figure
subplot(3,2,1),plot(y);,title('(a)原始平稳信号')
subplot(3,2,3),mesh(abs(GDfy)); xlabel('b'),ylabel('w'),zlabel('|Gabor(a,w)|');;title('(b)Gabor变换的频谱立体图')
subplot(3,2,4),imshow(abs(GDfy));;title('(c)Gabor变换的频谱频谱图图像')
subplot(3,2,5),mesh(abs(GDfy_shift)); xlabel('b'),ylabel('w'),zlabel('|Gabor(a,w)|');title('(d)对称平移的Gabor变换的频谱立体图')
subplot(3,2,6),imshow(abs(GDfy_shift));;title('(e)对称平移的Gabor变换的频谱立体图')
子程序
%输入变量a 为窗口控制参数,b为平移参数,x为函数变量
%输出变量gy为fun函数与窗口函数g的卷积
function gy = Gb_fun(a,b,x)
fun = @(x)0.*(x<0 |x>5)+x.*(x>=0 & x<1)+1.*(x>=1 & x<4) +(5-x).*(x>=4 & x<5);%分段函数
g = @(a,x)1./(2*(pi.*a).^(1/2)).*exp(-x.^2./(4*a));%窗口函数
gy =g(a,x-b.*ones(size(x))).*fun(x);% 窗口函数与原始信号的卷积