使用MATLAB直观地了解傅里叶变换
侧边栏壁纸
  • 累计撰写 18 篇文章
  • 累计收到 0 条评论

使用MATLAB直观地了解傅里叶变换

lycraft
2020-11-16 / 0 评论 / 37 阅读 / 正在检测是否收录...

使用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中是很难看出原本的波形成分的,但是经过傅里叶变换后的图形,我们就很容易看出原本的波形频率信息。

所以在一些场景下,频域比时域更适合进行研究。

0

评论

博主关闭了所有页面的评论