使用MATLAB直观地告诉你什么是傅里叶变换
要了解傅里叶变换 我们先来看看它在生活中都有哪些实际应用。
首先我们要知道 傅氏变换是将研究对象从时域变换到频域的
时域就是这样的一个东西啦。横坐标是t时间,经过傅里叶变换以后呢,横坐标就变成f频率了。在有些场景下,时域是很难分析出我们需要的结果的,但是从频域上,就可以分析出来。就比如声音信号和图像信号。
首先看声音的频谱吧
声音频谱的例子我准备用软件Au来演示,它里边有个功能就是FFT(快速傅里叶变换)
如图,声音的频谱中 x轴上离0点近的地方算作低频成分,随着f的增大,算是高频成分。我们知道 男生说话的声音频率比较低 所以如果给个高通滤波器,滤掉低频的男生声音部分,然后再傅里叶反变换回来,那么男生说话的声音就会被滤去。如果把低频成分向高频迁移的话,男生的声音就会变声成像女生的声音
然后看图像的频谱。
美颜就用到了傅里叶变换,拿磨皮来说。照片中 人脸上的痘痘其实就可以算是傅里叶变换后的高频部分了。我们滤去高频成分,剩下的就是人脸轮廓了。
可以看到,通过低通滤波器滤去高频之后 图像中的细节就没有了
上边三个图的MATLAB代码,摘自CSDN。
% 此处代码改自https://blog.csdn.net/u012900447/article/details/53126880
clc;clear;close all;
%高斯低通频域滤波
I=imread('http://www.lycraft.top/img/assassin.jpg');
subplot(1,2,1),imshow(I),title('原始图像');
I=double(I);
S=fftshift(fft2(I));
[M,N]=size(S);
n=2;
d0=30; %GLPF滤波,d0=5,15,30(程序中以d0=30为例)
n1=floor(M/2);
n2=floor(N/2);
for i=1:M
for j=1:N
d=sqrt((i-n1)^2+(j-n2)^2);
h=1*exp(-1/2*(d^2/d0^2));
S(i,j)=h*S(i,j);
end
end
S=ifftshift(S);
S=uint8(real(ifft2(S)));
subplot(1,2,2),imshow(S),title('高斯低通滤波图像');
%巴特沃斯低通频域滤波
I=imread('http://www.lycraft.top/img/assassin.jpg');
figure(2)
subplot(1,2,1),imshow(I),title('原始图像');
F=double(I); % 数据类型转换,MATLAB不支持图像的无符号整型的计算
G=fft2(F); % 傅立叶变换
G=fftshift(G); % 转换数据矩阵
[M,N]=size(G);
nn=2; % 二阶巴特沃斯(Butterworth)低通滤波器
d0=30; %截止频率为30
m=fix(M/2); n=fix(N/2);
for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
h=1/(1+0.414*(d/d0)^(2*nn)); % 计算低通滤波器传递函数
result(i,j)=h*G(i,j);
end
end
result=ifftshift(result);
Y2=ifft2(result);
Y3=uint8(real(Y2));
subplot(1,2,2),imshow(Y3),title('巴特沃斯低通滤波')
%巴特沃斯高通频域滤波
I=imread('http://www.lycraft.top/img/assassin.jpg');
figure(3)
subplot(1,2,1),imshow(I); title('原始图像');
F=double(I);% 数据类型转换,MATLAB不支持图像的无符号整型的计算
G=fft2(F);% 傅立叶变换
G=fftshift(G);%转换数据矩阵
[M,N]=size(G);
nn=2;% 二阶巴特沃斯(Butterworth)高通滤波器
d0=30;
m=fix(M/2);n=fix(N/2);
for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
if (d==0)
h=0;
else
h=1/(1+0.414*(d0/d)^(2*nn));% 计算传递函数
end
result(i,j)=h*G(i,j);
end
end
result=ifftshift(result);
J2=ifft2(result);
J3=uint8(real(J2));
subplot(1,2,2);imshow(J3); title('巴特沃斯高通滤波后图像'); % 滤波后图像显示
现在我们把两个正弦波合成以后,再把它进行傅里叶变换看看结果!
现在搞清楚了,感觉还是使用图像的方式最容易理解
先贴上MATLAB的代码
%% 清除缓存
clc;close all;clear;
%% init
Fs = 1000; %采样频率
t = 0:1/Fs:1; %从0时刻到1时刻采样1000次
L = length(t); %信号长度
%% 建立正弦信号
s1 = sin(2*pi*50*t);
subplot(5,1,1);plot(100*t(1:100),s1(1:100));
title('s1 = sin(2*pi*50*t)');xlabel('t(ms)');ylabel('s1');
s2 = sin(2*pi*150*t);
subplot(5,1,2);plot(100*t(1:100),s2(1:100));
title('s2 = sin(2*pi*150*t)');xlabel('t(ms)');ylabel('s2');
s12 = s1 + s2;
subplot(5,1,3);plot(100*t(1:100),s12(1:100));
title('s1+s2');xlabel('t(ms)');ylabel('s12');
%% 加入高斯噪声
s123 = s12 + randn(size(t));
subplot(5,1,4);plot(100*t(1:100),s123(1:100));
title('加入高斯白噪声后的两正弦函数相加');
%% 进行傅里叶变换
S1 = fft(s123);
n = 2^nextpow2(L);
S2 = abs(S1);
f = 0:1:(L-1);
subplot(5,1,5);plot(f(1:300),S2(1:300));
title('将带高斯噪声的两正弦函数相加的波形进行傅里叶变换')
xlabel('f(hz)')
这段代码就是两个正弦波s1和s2.(s1的频率为50hz,s2的频率为150hz。)s12就是s1和s2波形的相加。然后再在s12波形中加入高斯白噪声形成波形s123。
然后对带有高斯白噪声的波形s123进行傅里叶变换最后得到波形S2.这是可以看到经过傅里叶变换后的图形上边的50hz和150hz就是s1和s2的频率了。
我们在s123中是很难看出原本的波形成分的,但是经过傅里叶变换后的图形,我们就很容易看出原本的波形频率信息。
所以在一些场景下,频域比时域更适合进行研究。
评论