精华内容
下载资源
问答
  • 最早的电话使用的模拟信号原理图(1)...局限性就是当距离比较远 的时候就很难还原(电流的衰减大)傅里叶变换反傅里叶变换就解决了这个问题--使得数字信号通信代替模拟信号通信称为可能(1)任意一个是信号f(...

    最早的电话使用的模拟信号原理图


    (1)声音通过金属振动膜感应声波来影响磁场和电流,并将这种带有金属振动膜振动的“信息”的电流传递给另一端

    (2)另一端则进行反向工作,把不断变化的电流转化为电线圈中磁场的变化,使得金属振动膜发出同样的振动。

    局限性就是当距离比较远 的时候就很难还原(电流的衰减大)

    傅里叶变换和反傅里叶变换就解决了这个问题--使得数字信号通信代替模拟信号通信称为可能

    (1)任意一个是与信号f(t)经过傅里叶变换,都可以转化为多个余弦波叠加的形式

    一旦这个环节被我们掌控的话,语音信号就可以通过调制,由高频载波用极高的速率发送电信号和光信号。

    而在接收端,用滤波器过滤余弦波信号后,使用傅里叶反变换进行是与信号的还原


    对傅里叶的理解资源:

    深入浅出的讲解傅里叶变换(真正的通俗易懂) - CSDN博客
    https://blog.csdn.net/l494926429/article/details/51818012

    实验的来源:

    实验六傅里叶变换及其反变换_百度文库
    https://wenku.baidu.com/view/1ba79654178884868762caaedd3383c4bb4cb4bd.html


    代码:

    syms t v w x phase im re ;          %定义符号变量  
    a = input('请输入a=')
    f = exp(-a*abs(t));  %f(t) = exp(-2*t)*u(t)
    Fw = fourier(f);                     %求傅里叶变换
    subplot(311);
    ezplot(f);                          %绘制f(t)的时域波形
    axis([-1 2.5 0 1.1]);
    subplot(312);
    ezplot(abs(Fw));                    %绘制幅度谱
    im = imag(Fw);                      %计算F(w)的虚部
    re = real(Fw);                      %计算F(w)的实部
    phase = atan(im/re);                %计算相位谱
    subplot(313);
    ezplot(phase);                      %绘制相位谱


    syms t v w x phase im re ;          %定义符号变量  
    F(w) = 1/(1+w^2);                  %待求解变换的公式
    f = fourier(Fw,t);                 
    subplot(311);
    ezplot(f);                          %绘制f(t)的时域波形
    axis([-1 2.5 0 1.1]);
    subplot(312);
    ezplot(abs(Fw));                    %绘制幅度谱
    im = imag(Fw);                      %计算F(w)的虚部
    re = real(Fw);                      %计算F(w)的实部
    phase = atan(im/re);                %计算相位谱
    subplot(313);
    ezplot(phase);                      %绘制相位谱
    

    展开全文
  • 离散傅里叶变换与相关性计算

    千次阅读 2019-11-03 11:41:37
    文章目录离散傅里叶变换(DFT)相关性计算基本定义时域的卷积对应频域的相乘循环平移的傅里叶变换相关卷积通过傅里叶变换计算相关 本文主要介绍通过离散傅里叶变换的方法计算两个时间序列的相关性,不仅给出...

    离散傅里叶变换(DFT)与相关性计算


    本文主要介绍通过离散傅里叶变换的方法计算两个时间序列的相关性,不仅给出数学公式的推导,同时给出对应的 Matlab程序的验证。文中涉及到相关,卷积等基本概念,本文假设读者已经具备相关的数学基础。

    更改说明:由于CSDN的MD语法在解析部分公式时出现错误,所以讲部分公式改为图片。
    补充:关于循环卷积与频域的相乘关系

    基本定义

    首先,我们给出离散傅里叶变换,卷积与相关的定义:

    离散傅里叶变换

    长度为N的离散时间序列 x ( n ) x(n) x(n)的离散傅里叶变换的定义为:
    X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π n k / N X(k) = \sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N} X(k)=n=0N1x(n)ej2πnk/N
    对应的Matlab函数为fft(x)

    卷积

    长度为N的离散时间序列 x ( n ) x(n) x(n) y ( n ) y(n) y(n)的卷积定义为:
    z ( n ) = ∑ m = 0 N − 1 x ( m ) y ( n − m ) z(n) = \sum^{N-1}_{m=0} x(m)y(n-m) z(n)=m=0N1x(m)y(nm)
    得到的卷积 z ( n ) z(n) z(n)的长度为 2 N − 1 2N-1 2N1,n的取值范围是0到2N-2,卷积对应的Matlab函数为conv(x, y)。举个简单的例子:

    x = [1, 2, 3]
    y = [1, -1]
    z = conv(x, y) %  [1, 1, 1, -3]
    

    相关

    长度为N的离散时间序列 x ( n ) x(n) x(n)的相关定义为:
    z ( n ) = ∑ m = 0 N − 1 x ( m ) y ( m − n ) z(n) = \sum^{N-1}_{m=0} x(m)y(m-n) z(n)=m=0N1x(m)y(mn)
    同样,得到的相关 z ( n ) z(n) z(n)的长度为 2 N − 1 2N-1 2N1,卷积对应的Matlab函数为xcorr(x, y), 注意n的范围是 − ( N − 1 ) -(N-1) (N1) ( N − 1 ) (N-1) (N1)。举个简单的例子:

    x = [1, 2, 3]
    y = [1, -1, 1]
    z = xcorr(x, y) %  [1, 1, 2, -1, 3],对应的n是[-2, -1, 0, 1, 2]
    

    注意xcorr中要求两个时间序列长度相等,如果不相等,会自动将较短的时间序列补0。

    时域的卷积对应频域的相乘

    z ( n ) z(n) z(n)是时间序列 x ( n ) x(n) x(n) y ( n ) y(n) y(n)的卷积,长度为 2 N − 1 2N-1 2N1,则对应的离散傅里叶变换为:

    Z ( k ) = ∑ n = 0 2 N − 2 z ( n ) e − j 2 π n k / 2 N − 1 = ∑ n = 0 2 N − 2 ∑ m = 0 N − 1 x ( m ) y ( n − m ) e − j 2 π n k / 2 N − 1 \begin{aligned} Z(k) & = \sum^{2N-2}_{n=0} z(n)e^{-j2\pi nk/{2N-1}} \\ & = \sum^{2N-2}_{n=0} \sum^{N-1}_{m=0} x(m)y(n-m)e^{-j2\pi nk/{2N-1}} \end{aligned} Z(k)=n=02N2z(n)ej2πnk/2N1=n=02N2m=0N1x(m)y(nm)ej2πnk/2N1
    交换求和顺序得
    Z ( k ) = ∑ m = 0 N − 1 x ( m ) ∑ n = 0 2 N − 2 y ( n − m ) e − j 2 π n k / 2 N − 1 = ∑ m = 0 N − 1 x ( m ) e − j 2 π m k / 2 N − 1 ∑ n = 0 2 N − 2 y ( n − m ) e − j 2 π ( n − m ) k / 2 N − 1 = ∑ m = 0 N − 1 x ( m ) e − j 2 π m k / 2 N − 1 ∑ n = 0 N − 1 y ( n ) e − j 2 π n k / 2 N − 1 \begin{aligned} Z(k) & = \sum^{N-1}_{m=0} x(m) \sum^{2N-2}_{n=0} y(n-m)e^{-j2\pi nk/{2N-1}} \\ & = \sum^{N-1}_{m=0} x(m)e^{-j2\pi mk/{2N-1}} \sum^{2N-2}_{n=0} y(n-m)e^{-j2\pi (n-m)k/{2N-1}} \\ & = \sum^{N-1}_{m=0} x(m)e^{-j2\pi mk/{2N-1}} \sum^{N-1}_{n=0} y(n)e^{-j2\pi nk/{2N-1}} \end{aligned} Z(k)=m=0N1x(m)n=02N2y(nm)ej2πnk/2N1=m=0N1x(m)ej2πmk/2N1n=02N2y(nm)ej2π(nm)k/2N1=m=0N1x(m)ej2πmk/2N1n=0N1y(n)ej2πnk/2N1
    注意等式(7)中第二项的变量替换。由于 Z ( k ) Z(k) Z(k)序列的长度为 2 N − 1 2N-1 2N1,而 X ( k ) X(k) X(k) Y ( k ) Y(k) Y(k)的长度都为 N N N,为了长度上一致,我们将 x ( n ) x(n) x(n) y ( n ) y(n) y(n)分别补0,使得最终的序列长度为 2 N − 1 2N -1 2N1,分别记作 x p ( n ) x_p(n) xp(n) y p ( n ) y_p(n) yp(n),以 x p ( n ) x_p(n) xp(n)为例对应的数学描述为:

    x p ( n ) = { x ( n ) n ≤ N − 1 0 n > N − 1 x_p(n) = \left\{ \begin{array}{rl} x(n) & n \leq N-1\\ 0 & n>N-1 \end{array} \right. xp(n)={x(n)0nN1n>N1
    则其对应的傅里叶变换为:
    X p ( k ) = ∑ n = 0 N − 1 x p ( n ) e − j 2 π n k / 2 N − 1 X_p(k) = \sum_{n=0}^{N-1}x_p(n)e^{-j2\pi nk/2N-1} Xp(k)=n=0N1xp(n)ej2πnk/2N1
    联系等式(8)可以得到
    Z ( k ) = X p ( k ) Y p ( k ) Z(k) = X_p(k)Y_p(k) Z(k)=Xp(k)Yp(k)
    我们使用一个简单的Matlab程序进行验证。

    N = 20;
    x = randn(1, N);
    y = randn(1, N);
    
    z = conv(x, y);
    
    xp = [x, zeros(1, N-1)]; %补0
    yp = [y, zeros(1, N-1)]; %补0
    
    zk = fft(z);
    
    xpk = fft(xp);
    ypk = fft(yp);
    zke = xpk.*ypk;%频域相乘
    
    real_err = max(abs(real(zk)-real(zke)));
    imag_err = max(abs(imag(zk) - imag(zke)));
    
    fprintf('实补误差 %e\n虚部误差 %e\n',real_err, imag_err);
    
    %% 程序输出结果
    % 实部误差 1.421085e-14
    % 虚部误差 8.881784e-15
    

    其他的文献补0的个数都选择N而非N-1,背后的道理是一样的。选择N的好处是 X p ( k ) X_p(k) Xp(k) X ( k ) X(k) X(k)有一个一致的对应关系,即 X p ( 2 k ) = X ( k ) X_p(2k) = X(k) Xp(2k)=X(k)

    上面的定理可以简单描述为时域的卷积对应频域的乘积,一个典型应用加速卷积的计算。我们知道,如果直接套用卷积的定义,那么时间复杂度是 O ( n 2 ) O(n^2) O(n2),而快速傅里叶变换的时间复杂度是 O ( n log ⁡ n ) O(n\log n) O(nlogn),一个直接的思路就是利用傅里叶变换计算卷积,即首先根据 X p ( k ) X_p(k) Xp(k) Y p ( k ) Y_p(k) Yp(k)得到时间 Z ( k ) Z(k) Z(k),然后进行傅里叶反变换得到卷积 z ( n ) z(n) z(n)

    循环卷积与相乘的关系

    在图像处理中,图像与卷积核的运算通常需要保持图像大小不变,这种情况下时域的卷积与频域的乘积的对应关系是否仍然保持呢?
    Matlab对应的卷积函数为imfilter(u, v,'circular', 'conv'),给出简单的Matlab函数进行验证。

    a=[1,2,3,4,5,6,7];
    b = [1,2,3];
    d = imfilter(a,b,'circular','conv');
    % 结果 d => [ 11     8    12    16    20    24    21]
    

    了解了循环卷积的运算过程,我们进行公式的推导。首先,对循环卷积进行表示。对于长度为N的序列 x ( n ) x(n) x(n)其中,n=0, 1, 2, …N-1。以及长度为M的序列 h ( n ) h(n) h(n),这里对于 h ( n ) h(n) h(n),如果长度为奇数2M+1,则对应的n的取值范围为-M, -(M-1), … 0, …, (M-1), M。如果长度为偶数2M,则对应的n的范围表示为-M, …M-1。则循环卷积y(n)可以表示为:
    y ( n ) = ∑ m = 0 N − 1 x ( m ) h ) ( n − m ) , n = 0 , 1 , 2 , . . . N − 1 y(n) = \sum_{m=0}^{N-1}x(m)h)(n-m), n=0,1,2,...N-1 y(n)=m=0N1x(m)h)(nm),n=0,1,2,...N1
    根据上述公式,进行形同推导,则可以得到:
    Y ( k ) = ∑ m = 0 N − 1 x ( m ) ∑ n = 0 N − 1 h ( n − m ) e − j 2 π n k / N = ∑ m = 0 N − 1 x ( m ) e − j 2 π m k / N ∑ n = 0 N − 1 h ( n − m ) e − j 2 π ( n − m ) k / N \begin{aligned} Y(k) & = \sum^{N-1}_{m=0} x(m) \sum^{N-1}_{n=0} h(n-m)e^{-j2\pi nk/N} \\ & = \sum^{N-1}_{m=0} x(m)e^{-j2\pi mk/N} \sum^{N-1}_{n=0} h(n-m)e^{-j2\pi (n-m)k/N} \\ \end{aligned} Y(k)=m=0N1x(m)n=0N1h(nm)ej2πnk/N=m=0N1x(m)ej2πmk/Nn=0N1h(nm)ej2π(nm)k/N
    定义 h p ( n ) h_p(n) hp(n)为将 h ( n ) h(n) h(n)补0到长度为N,然后将向左循环平移M位。即:
    h p ( n ) = { h ( n ) 0 ≤ n ≤ M − 1 x ( n − m + N ) n ≥ n − M 0 o t h e r w i s e h_p(n) = \left\{ \begin{array}{rl} h(n) & 0\leq n \leq M-1\\ x(n-m+N) & n \geq n-M\\ 0 & otherwise \end{array} \right. hp(n)=h(n)x(nm+N)00nM1nnMotherwise
    Y ( k ) = X ( k ) H p ( k ) Y(k)=X(k)H_p(k) Y(k)=X(k)Hp(k),我们通过一个Matlab程序进行验证。

    clc
    clear
    x=[1,2,3,4,5,6,7];
    h = [1,2,-1];
    
    y = imfilter(x, h, 'circular', 'conv'); %-3     6     8    10    12    14     9
    
    
    hp = [h, zeros(1, length(x)-length(h))]; %补0
    m = floor(length(h)/2);
    hp = circshift(hp, -m); %向左循环平移
    Yk = fft(x).*fft(hp);
    
    ye = ifft(Yk);  %  -3.0000    6.0000    8.0000   10.0000   12.0000   14.0000    9.0000
    % 计算结果可知,y与ye的结果一致(忽略舍入误差)
    

    循环平移的傅里叶变换

    循环平移的本质是将序列 x ( n ) x(n) x(n)看做周期为N的周期序列,也就是N到2N-1的取值与0到N-1的取值有一一对应的关系。这里我们给出证明:

    将序列 x ( n ) x(n) x(n)向右做循环平移m位,得到y(n),则:
    y ( n ) = { x ( n − m ) n ≥ m x ( n − m + N ) n < m y(n) = \left\{ \begin{array}{rl} x(n-m) & n \geq m\\ x(n-m+N) & n < m \end{array} \right. y(n)={x(nm)x(nm+N)nmn<m
    则对应的傅里叶变换为:
    Y ( k ) = ∑ n = 0 N − 1 y ( n ) e − j 2 π n k / N = ∑ n = m N − 1 x ( n − m ) e − j 2 π n k / N + ∑ n = 0 m − 1 x ( n − m + N ) e − j 2 π n k / N = ∑ n = m N − 1 x ( n − m ) e − j 2 π ( n − m ) k / N e − j 2 π m k / N + ∑ n = 0 m − 1 x ( n − m + N ) e − j 2 π ( n − m + N ) k / N e − j 2 π ( m − N ) k / N = ∑ n = 0 N − 1 x ( n ) e − j 2 π n k / N e − j 2 π m k / N = X ( k ) e − j 2 π m k / N \begin{aligned} Y(k) &= \sum_{n=0}^{N-1} y(n)e^{-j2\pi nk/N} \\ &= \sum_{n=m}^{N-1} x(n-m)e^{-j2\pi nk/N} + \sum_{n=0}^{m-1} x(n-m+N)e^{-j2\pi nk/N} \\ &= \sum_{n=m}^{N-1} x(n-m)e^{-j2\pi (n-m)k/N}e^{-j2\pi mk/N} + \sum_{n=0}^{m-1} x(n-m+N)e^{-j2\pi (n-m+N)k/N}e^{-j2\pi (m-N)k/N} \\ &= \sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N}e^{-j2\pi mk/N} \\ &= X(k)e^{-j2\pi mk/N} \end{aligned} Y(k)=n=0N1y(n)ej2πnk/N=n=mN1x(nm)ej2πnk/N+n=0m1x(nm+N)ej2πnk/N=n=mN1x(nm)ej2π(nm)k/Nej2πmk/N+n=0m1x(nm+N)ej2π(nm+N)k/Nej2π(mN)k/N=n=0N1x(n)ej2πnk/Nej2πmk/N=X(k)ej2πmk/N

    相关与卷积

    熟悉卷积的读者应该知道卷积的计算过程分为两步:翻转和平移,与相关的计算相比,只是多了一步翻转操作。如果直接序列进行翻转,是否能用卷积代替相关操作呢?答案是是显然的,我们用一个Matlab程序做证明。

    x = [1,2,1];
    y = [1,2,3];
    
    z1 = xcorr(x, y); %相关
    z2 = conv(x, flip(y)); %先翻转,然后卷积
    
    z1 % => [3, 8, 8, 4, 1]
    z2 % => [3, 8, 8, 4, 1]
    

    这里不做数学的证明,有兴趣的读者可以自己证明。证明的过程中注意相关计算结果的第一项对应的n是 − ( N − 1 ) -(N-1) (N1),而卷积操作的第一项对应的n是0。

    通过傅里叶变换计算相关

    如果直接利用相关的公式,时间复杂度是 O ( n 2 ) O(n^2) O(n2),而快速傅里叶变换的时间复杂度是 O ( n log ⁡ n ) O(n\log n) O(nlogn),那么能否使用fft方法来计算相关,从而提高运算速度呢?上一节我们已经知道相关与翻转后做卷积是等价的,那么一个直接的思路是:

    1. 将序列 y ( n ) y(n) y(n)翻转,得到 u ( n ) u(n) u(n),然后补0,得到 u p ( n ) u_p(n) up(n),进行傅里叶变换得到 U p ( k ) U_p(k) Up(k)
    2. 将序列 x ( n ) x(n) x(n)补0,得到 x p ( n ) x_p(n) xp(n),进行傅里叶变换得到 X p ( k ) X_p(k) Xp(k)
    3. 计算 Z ( k ) = X ( p ) U p ( k ) Z(k) = X_(p)U_p(k) Z(k)=X(p)Up(k),通过傅里叶反变换得到 z ( n ) z(n) z(n)

    上述计算的缺点是多了一个中间变量 u ( n ) u(n) u(n),那么能否用 Y p ( k ) Y_p(k) Yp(k) 来计算 U p ( k ) U_p(k) Up(k)呢?我们进行一个简单的推导。
    U p ( k ) = ∑ n = 0 N − 1 u ( n ) e − j 2 π n k / 2 N − 1 = ∑ n = 0 N − 1 y ( N − 1 − n ) e − j 2 π n k / 2 N − 1 = ∑ n = 0 N − 1 y ( N − 1 − n ) e − j 2 π ( n + 1 − N ) k / 2 N − 1 e − j 2 π ( N − 1 ) k / 2 N − 1 = { ∑ n = 0 N − 1 y ( N − 1 − n ) e − j 2 π ( N − 1 − n ) k / 2 N − 1 } ∗ e − j 2 π ( N − 1 ) k / 2 N − 1 = Y p ∗ ( k ) e − j 2 π ( N − 1 ) k / 2 N − 1 \begin{aligned} U_p(k) & = \sum_{n=0}^{N-1}u(n)e^{-j2\pi nk/2N-1} \\ & = \sum_{n=0}^{N-1}y(N-1-n)e^{-j2\pi nk/2N-1} \\ & = \sum_{n=0}^{N-1}y(N-1-n)e^{-j2\pi (n+1-N)k/2N-1}e^{-j2\pi (N-1)k/2N-1} \\ & = \left\{\sum_{n=0}^{N-1}y(N-1-n)e^{-j2\pi (N-1-n)k/2N-1}\right\}^*e^{-j2\pi (N-1)k/2N-1} \\ & = Y^*_p(k)e^{-j2\pi (N-1)k/2N-1} \end{aligned} Up(k)=n=0N1u(n)ej2πnk/2N1=n=0N1y(N1n)ej2πnk/2N1=n=0N1y(N1n)ej2π(n+1N)k/2N1ej2π(N1)k/2N1={n=0N1y(N1n)ej2π(N1n)k/2N1}ej2π(N1)k/2N1=Yp(k)ej2π(N1)k/2N1
    其中 ∗ * 表示共轭运算,对应的matlab函数为conj(X)

    也就是说如果令 Z 1 ( k ) = X p ( k ) Y p ∗ ( k ) Z_1(k) = X_p(k)Y_p^*(k) Z1(k)=Xp(k)Yp(k),对应的傅里叶反变换为 z 1 ( n ) z_1(n) z1(n),则 Z ( k ) = Z 1 ( k ) e − j 2 π ( N − 1 ) k / 2 N − 1 Z(k) = Z_1(k)e^{-j2\pi (N-1)k/2N-1} Z(k)=Z1(k)ej2π(N1)k/2N1,根据循环平移的傅里叶变换易知 z ( n ) = z 1 ( n − ( N − 1 ) ) z(n) = z1(n-(N-1)) z(n)=z1(n(N1))

    因此,通过傅里叶变换计算相关运算的思路为:

    1. x ( n ) x(n) x(n) y ( n ) y(n) y(n)分别补N-1个0,得到 x p ( n ) x_p(n) xp(n) y p ( n ) y_p(n) yp(n),进行快速傅里叶变换得到 X p ( k ) X_p(k) Xp(k) Y p ( k ) Y_p(k) Yp(k)
    2. 计算 Z 1 ( k ) = X ( k ) Y ∗ ( k ) Z_1(k) = X(k)Y^*(k) Z1(k)=X(k)Y(k),并进行傅里叶反变换得到 z 1 ( n ) z_1(n) z1(n)
    3. 进行循环平移得到 z ( n ) z(n) z(n)

    同之前讨论过的一样,补0的个数只要大于等于N-1不影响计算结果,很多资料中都补N个0.

    最后我们通过一个Matlab程序对比一下4种计算相关性的方法。

    % 通过4种方式计算相关
    %% 数据生成
    N = 30;
    x = randn(1, N);
    y = [randn(1, 4), x(1:N-4)];
    %y = randn(1, N);
    
    %% 方法1, 通过Matlab自带函数xcorr进行计算
    
    z1 = xcorr(x, y);
    
    %% 方法2,通过翻转,求卷积
    
    z2 = conv(x, flip(y));
    
    %% 方法3,通过翻转,傅里叶变换计算
    u = flip(y);
    xp = [x, zeros(1, N-1)];
    up = [u, zeros(1, N-1)];
    
    zk = fft(xp).*fft(up);
    z3 = ifft(zk);
    
    %% 方法4,傅里叶变换方法
    xp = [x, zeros(1, N-1)];
    
    yp = [y, zeros(1, N-1)];
    
    z1k = fft(xp).*conj(fft(yp));
    z1n =  ifft(z1k);
    z4 = circshift(z1n, N-1); %需要平移操作
    
    %% 计算结果的可视化
    lags = -(N-1):N-1;
    
    
    scalediff = 5;
    plot(lags, z1);
    hold on;
    plot(lags, z2 + scalediff*1);
    plot(lags, z3 + scalediff*2);
    plot(lags, z4 + scalediff*3);
    
    plot([-4, -4], ylim, 'r:','LineWidth', 1.5); % 平移量
    hold off
    legend('方法一','方法二','方法三', '方法四');
    
    

    结果:
    在这里插入图片描述

    展开全文
  • 计算信号的傅里叶变换,信号的反傅里叶变换,周期信号的展开,
  • 3.1 Python图像的频域图像增强-图像的傅里叶变换反变换 文章目录3.1 Python图像的频域图像增强-图像的傅里叶变换反变换1 算法原理2 代码3 效果 1 算法原理 图像的傅里叶变换反变换(需要考虑图像旋转、平移时...

    3.1 Python图像的频域图像增强-图像的傅里叶变换和反变换

    1 算法原理

    图像的傅里叶变换和反变换(需要考虑图像旋转、平移时的变换)

    图像的傅里叶变换和反变换

    二维离散傅里叶变换(Two-Dimensional Discrete Fourier Transform)常用于图像处理中,对图像进行傅里叶变换后得到其频谱图。频谱图中频率高低表征图像中灰度变化的剧烈程度。图像中边缘和噪声往往是高频信号,而图像背景往往是低频信号。我们在频率域内可以很方便地对图像的高频或低频信息进行操作,完成图像去噪,图像增强,图像边缘提取等操作。图像(M*N)的二维离散傅立叶变换可以将图像由空间域变换到频域中去,空间域中用 x,y 来表示空间坐标,频域由 u,v 来表示频率,二维离散傅立叶变换的公式如下:

    image-20210710092425135

    图像长 M,高 N。F(u,v)表示频域图像,f(x,y)表示时域图像。u 的范围为[0,M-1],v 的范围为[0,N-1]。

    对二维图像进行傅里叶逆变换式子如下:

    image-20210710092439193

    图像长 M,高 N。f(x,y)表示时域图像, F(u,v)表示频域图像。x 的范围为[0,M-1],y 的范围为[0,N-1]。

    在 python 中,numpy 库的 fft 模块有实现好了的二维离散傅立叶变换函数,函数是 fft2,输入一张灰度图,输出经过二维离散傅立叶变换后的结果,但是具体实现并不是直接用上述公式,而是用快速傅立叶变换。结果需要通过使用 abs 求绝对值才可以进行可视化,但是视觉效果并不理想,因为傅立叶频谱范围很大,所以要用 log 对数变换来改善视觉效果。

    在使用 log 函数的时候,要写成 log(1 + x) 而不是直接用 log(x),这是为了避开对 0 做对数处理。另外,图像变换的原点需要移动到频域矩形的中心,所以要对 fft2 的结果使用 fftshift 函数。最后也可以使用 log 来改善可视化效果。

    2 代码

    运行代码说明

    1.要改变代码中的图片地址(地址不能有中文)

    更改put(path)函数中的路径put(r'../image/image1.jpg')

    2.注意最后的plt.savefig('1.new.jpg')是保存plt图像,如果不使用可以注释掉

    import os
    import numpy as np
    import cv2
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    
    def put(path):
        img = cv2.imread(path, 1)
        # img = cv2.imread(os.path.join(base, path), 1)
        img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        rows, cols = img.shape[:2]
        # 傅里叶变换
        dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
        # 将频谱低频从左上角移动至中心位置
        dft_shift = np.fft.fftshift(dft)
        # 频谱图像双通道复数转换为0-255区间
        res1 = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
        # 傅里叶逆变换
        ishift = np.fft.ifftshift(dft_shift)
        iimg = cv2.idft(ishift)
        res2 = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])
        # 图像顺时针旋转60度
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), -60, 1)
        rot = cv2.warpAffine(img, M, (rows, cols))
        # 傅里叶变换
        dft3 = cv2.dft(np.float32(rot), flags=cv2.DFT_COMPLEX_OUTPUT)
        dft_shift3 = np.fft.fftshift(dft3)
        res3 = 20 * np.log(cv2.magnitude(dft_shift3[:, :, 0], dft_shift3[:, :, 1]))
        # 傅里叶逆变换
        ishift3 = np.fft.ifftshift(dft_shift3)
        iimg3 = cv2.idft(ishift3)
        res4 = cv2.magnitude(iimg3[:, :, 0], iimg3[:, :, 1])
    
        # 图像向右平移
        H = np.float32([[1, 0, 200], [0, 1, 0]])
        tra = cv2.warpAffine(img, H, (rows, cols))
        # 傅里叶变换
        dft2 = cv2.dft(np.float32(tra), flags=cv2.DFT_COMPLEX_OUTPUT)
        dft_shift2 = np.fft.fftshift(dft2)
        res5 = 20 * np.log(cv2.magnitude(dft_shift2[:, :, 0], dft_shift2[:, :, 1]))
        # 傅里叶逆变换
        ishift2 = np.fft.ifftshift(dft_shift2)
        iimg2 = cv2.idft(ishift2)
        res6 = cv2.magnitude(iimg2[:, :, 0], iimg2[:, :, 1])
        # 输出结果
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.subplot(331), plt.imshow(img, plt.cm.gray), plt.title('原图灰度图像'),plt.axis('off')
        plt.subplot(332),plt.imshow(res1,plt.cm.gray),plt.title('傅里叶变换'),plt.axis('off')
        plt.subplot(333),plt.imshow(res2,plt.cm.gray),plt.title('傅里叶反变换'),plt.axis('off')
        plt.subplot(334),plt.imshow(rot,plt.cm.gray),plt.title('图像旋转'),plt.axis('off')
        plt.subplot(335),plt.imshow(res3,plt.cm.gray),plt.title('傅里叶变换'),plt.axis('off')
        plt.subplot(336),plt.imshow(res4,plt.cm.gray),plt.title('傅里叶反变换'),plt.axis('off')
        plt.subplot(337),plt.imshow(tra,plt.cm.gray),plt.title('图像平移'),plt.axis('off')
        plt.subplot(338),plt.imshow(res5,plt.cm.gray),plt.title('傅里叶变换'),plt.axis('off')
        plt.subplot(339),plt.imshow(res6,plt.cm.gray),plt.title('傅里叶反变换'),plt.axis('off')
    
        # plt.savefig('1.new.jpg')
        plt.show()
    
    # 处理函数,要传入路径
    put(r'../image/image3.jpg')
    

    3 效果

    image-20210710092845353

    展开全文
  • 计算傅里叶变换 傅里叶变换 鉴于这种想法,任何信号,当然任何周期性信号,都可以由一系列正弦曲线组成,我们将开始从级数(Series)的概念转向连续信号的概念。我们将要讨论一些东西,这些东西可以让我们知道图像...

    目录

    傅里叶变换

    计算傅里叶变换


    傅里叶变换

    鉴于这种想法,任何信号,当然任何周期性信号,都可以由一系列正弦曲线组成,我们将开始从级数(Series)的概念转向连续信号的概念。我们将要讨论一些东西,这些东西可以让我们知道图像中任意给定频率的功率有多大。也就是说,我们想把图像从某样东西上变换过来,该东西是时间的函数,或者仅仅是空间的函数,从而知道它的频率是多少。这种变换被称为什么?它叫做傅里叶变换(Fourier Transform)

    现在我们做了傅里叶级数(Fourier series),然后我们进入傅里叶变换(Fourier Transform),再然后我们将进行离散傅里叶变换(Discrete Fourier Transform, DFT),这将使我们得到那个图像的东西。

    所以我们想了解频率,频率通常被写成\omega(Omega),我们的信号。我们想要参数化,我们想要变换信号,而不是x,在空间中是通过\omega(Omega)。我们通过所谓的傅里叶变换(Fourier Transform)来实现这一点,傅里叶变换会拿一些空间信号,

    并且它将为我们提供一些特性,这些特性可以给每个给定\omega(Omega)的相位\varphi 和 幅度A进行编码。

    好的,所以我们需要给定每个 \omega (Omega)的相位\varphi 和 幅度A。特别是每个\omega (Omega)从0到无穷大,实际上它是 从负无穷大到无穷大(-\infty+\infty)。这个想法是f(\omega)保存了正弦曲线相应的幅度A 和 相位\varphi。因此,这个f 必须拥有幅度A 和 相位\varphi

    f 如何保留振幅和相位?复数的把戏(Complex number trick!)。注意,我没有说过一个单独的数字。因为请记住,f(x)只是F的值,所以它是一个数字,f(\omega)将是幅度和相位,我们该怎么做?非常好,实际上,大写字母F实际上是一个复数。现在,希望您能记住关于复数的一些知识。如果没有,我会给你一个非常简短的回顾。所以基本上, f(\omega)由两部分组成:一个实部(Real Part),另一个是虚部(Imaginary Part)

    所以请记住:A+Bi ,在这种情况下,  f(\omega)的实部加上虚部,也许你还记得复数的大小就是两个元素之和的平方根。通常是\sqrt{A^2+B^2},这里是 \sqrt{R^2+I^2}

    R代表实部,I是虚部,相位之间的关系将以这种方式编写,相位\varphi是 虚部 比 实部 的反正切。

    这里遗漏的点是:实部归因于 偶函数(Even),虚部归因于 奇函数(Odd),

    所有这意味着:如果我在这里有一个函数,是我的cosign函数,这是关于原点对称的。

    正弦曲线,我要把它们拧紧,因为它必须到这里,

    然后在那里,

    再然后,所以它必须是0(zero),在这里,这是我的正弦曲线,因为当那个东西高时,它必须降到0,奇怪的是,Sin (- x) 等于 -Sin (x) ,而 Cos (- x) 等于 Cos (x)。

    所以最后的这个部分(如图),

    实部是余弦部分(Cosin),

    虚部是正弦部分(Sine)。

    计算傅里叶变换

    当我们进行傅里叶变换时,我们所做的只是计算基集。我会告诉你我的意思。看到这个丑陋,丑陋的整体吗?

    好吧,我得到的是,我已经得到了一个函数的正弦值,另一个函数的正弦值。,

    两个不同的频率,a 和 b

    我将这一点从负无穷大加到正无穷大,

    我声称如果 a 不等于 b,它等于0。

    为什么这是正确的呢?有猜测吗? 让我们直观地思考一下,如果有一个正弦曲线,

    现在我必须做一个不同的频率,这很容易,

    天啊,真漂亮!哈哈? 在本质上,在某个点上,当这个是正的这个是相同的正值,我会试着找到,希望有一些是正的,

    我必须补上一个点。对,让我们假装红色的东西在这里下来就像那样 一点点。

    当红色的点是正数时,黄色的点是负数,

    所以这两个数的乘积会相互抵消,这是一种波浪的方式表示整个积分是0,只要 a \neq b

    但是你可能会问当a=b时会发生什么?

    我们先假设它们处于相同的相位,完全相同的相位。好吧,它就是Sin²。Sin²是正的,在无穷远处求和,得到什么?你得到无穷大。所以,如果你有两个频率相同的正弦波,你会得到一个无穷大的值。除非它们完全不相合,在这种情况下是0。换句话说,Sin乘以Cos继续下去。但最基本的思路是,如果我整合一个函数,如果它是由一个与\omega不同的正弦曲线组成的,当我采用\omega的正弦并且我完全间隔时,我什么也得不到。但是如果它是等于那个正弦,我会得到无穷大。而且这将是这样写的。

    所以让我们做一个简单的例子。假设有一个简单的函数:f(x)=cos(2 \pi \omega x)

    然后我们选一些频率,我们把它叫做u,简单点。

    如果我取这个积分,Okay。那么如果 u 等于那个相同的 \omega 那将是无穷的,否则它将是零。

    所以看起来像这样, Okay。我们只有这两个成员脉冲。这里只有这两个无穷大的尖峰。

    这被称为与余弦对应的脉冲。你可以看到它们是正的,这是因为它是余弦。

    如果我们有一个正弦,因为sin (- x) = - sin (x),它看起来就像在这里的向上,

    和另一个将会在这里向下。

    这就是正弦,Okay。这就是虚部部分。两个公式如下:

    我们只是在计算一个基集来说明这个正弦曲线有多少。我们可以为所有频率做到这一点。但是你可能会问,我们不是也必须为所有相位做到这一点吗?

    答案是否定的。最酷的事情之一是,如果我有cos(\omega x) 和 sin(\omega x)。所以1是90度 或 两相相变的功率。我可以通过它们的线性组合得到任意相位。你可以证明,这是真的,我们可以做一个小演示。但基本上,如果我有一个任意相位的正弦曲线,如果你告诉我它对余弦的积分是多少、它对正弦的积分是多少,我就能告诉你这个正弦曲线的相位是多少。


    复数的基础:

    https://blog.csdn.net/sw3300255/article/details/83149483


    ——学会编写自己的代码,才能练出真功夫。

    展开全文
  • python实现傅里叶变换反变换

    千次阅读 2019-04-16 00:05:02
    import numpy as np from math import * x=[1,2,3,4,5] ...'''傅里叶变换''' def fft(x): x = np.asarray(x, dtype=float) N = x.shape[0] n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * n...
  • 提出了分数傅里叶变换计算全息图(FRTCGH).采用分数傅里叶变换的快速算法和罗曼Ⅲ型迂回位相编码方法设计并...分析了FRTCGH和再现图像随分数阶变化的规律,讨论了分数傅里叶变换计算全息图菲涅耳计算全息图之间的关系.
  • 提出了双重分数傅里叶变换计算全息。在这种方法中,将两个图像的信息分别经不同阶的分数傅里叶变换后,记录在同一张分数傅里叶变换计算全息图上。它需要两个特定的分数傅里叶变换系统才能再现出所记录的图像信息,利用...
  • Fw=fourier(y):连续傅里叶变换函数 ft=ifourier(Fw,t):连续傅里叶变换函数 傅里叶变换举例如下: syms t; figure; y = (sin(2pi(t-1)))/(pi*(t-1)); ezplot(y);%创建了函数y并画图现实 Fw=fourier(y);%Fw即为...
  • 十分简明易懂的FFT(快速傅里叶变换

    万次阅读 多人点赞 2018-08-07 11:38:56
    快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶...
  • OpenCV中傅里叶变换反变换的运用

    千次阅读 2014-05-17 13:07:28
    //显示结果表明只取出反变换后的16%信息,在利用idft反变换,发现原图相差不太大 //这是图像处理中利用正交变换(包括傅里叶变换)的能量集中特性(图像有损压缩中有利用) Mat partFrequencyImg; ...
  • 初步了解快速傅里叶变换反变换的处理过程。 C++语言编写 #include <math.h> //编写头文件,还可采用#include "***"方式编写头文件 #include <malloc.h>//添加需要使用到的头文件 #define pi (double...
  • 傅里叶变换(二维离散傅里叶变换)

    万次阅读 多人点赞 2018-06-15 22:22:35
    离散二维傅里叶变换一常用性质: 可分离性、周期性和共轭对称性、平移性、旋转性质、卷积相关定理;(1)可分离性: 二维离散傅里叶变换DFT可分离性的基本思想是DFT可分离为两次一维DFT。因此可以用通过计算两次...
  • 一开始在思考怎么样进行傅里叶变换,要怎么进行复数的存储、计算、和表达呢,想了好久没能想出好的办法,百度了一下,没有找到好的文章来说明怎么存储,倒是找到了好多介绍怎么用opencv进行傅里叶变换的文章,详细一...
  • N点离散傅里叶变换同时计算两个N点实序列的离散傅里叶变换
  • 条纹投影轮廓测量中不连续点的相位计算傅里叶变换,开窗傅里叶变换和小波变换方法的比较
  • 1.图像频域处理的意义 在图像处理和分析中,经常会将图像从图像空间转换到其他空间中,并利用这些空间的特点进行对转换后图像进行分析处理,然后再将处理后的图像...傅里叶变换计算量较大,人们为了提高速度,
  • DFT的matlab源代码 #编制和运行离散傅里叶变换(DFT)的C语言程序,要求实现傅里叶正、反变换并获得正确的计算结果.
  • 傅里叶变换到加窗傅里叶变换

    万次阅读 多人点赞 2017-12-10 14:00:49
    傅里叶变换到加窗傅里叶变换我们都做了什么
  • OpenCV中的图像变换——傅里叶变换

    万次阅读 多人点赞 2021-07-22 20:11:24
    这篇博客将介绍OpenCV中的图像变换,包括用Numpy、OpenCV计算图像的傅里叶变换,以及傅里叶变换的一些应用;
  • 这是在一个论坛看的帖子, ...但若你问我,没学好傅里叶变换,能否操作(编程)小波变换,或是没学好第一代小波,能否操作二代小波变换,答案是肯定的。 一、基的概念 我们要明确的是基的概念。两者都是基,信
  • 傅里叶变换与图像处理

    千次阅读 2016-04-21 15:02:55
    本文主要介绍的是傅里叶变换在图像处理当中的应用。本文参考:http://blog.csdn.net/masibuaa/article/details/6316319 第一部分介绍傅里叶变换傅里叶变换是将时域信号分解为不同频率的正弦信号或余弦函数叠加之和...
  • 傅里叶变换与EEG傅里叶变换处理

    千次阅读 2019-07-31 20:31:25
    选取数据集Music BCI (006-2015)...由于数据量较大,每隔1000个点进行采样,对采样点进行傅里叶变换。 dataFir=load('musicbci_VPaak.mat');%读取数据 dataf=dataFir.data.X'; data1fft=dataf(1,:); [hei,wid]=s...
  • Q:简述计算机三大变换的联系和区别 (傅里叶变换 拉普拉斯变换 z变换)(1) 傅里叶变换定义:表示能将满足一定条件的某个函数表示成三角函数(正弦和/或余弦函数)或者它们的积分的线性组合。傅立叶变换是一种...
  • 傅里叶变换(非周期函数):傅里叶变换傅里叶级数变换而来,且傅里叶变换的应用不仅限于周期函数,也适用于非周期函数。 1. 三角函数系的正交性 定义三角函数系为这样一个集合:; 三角函数系的正交性是指:从...
  • 傅里叶变换&短时傅里叶变换&小波变换

    万次阅读 多人点赞 2016-05-03 11:26:36
    傅里叶变换&短时傅里叶变换&小波变换
  • 傅里叶变换 一维离散傅里叶变换

    万次阅读 热门讨论 2019-11-06 21:08:43
    DFT:(Discrete Fourier Transform)离散傅里叶变换傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域采样。在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列...
  • 实验六 傅里叶变换及其反变换6.1实验目的1.学会运用MATLAB 求连续时间信号的傅里叶变换;2.学会运用MATLAB 求连续时间信号的傅里叶反变换;3.学会运用MATLAB 求连续时间信号的频谱图。6.2实验原理及实例分析1...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 26,173
精华内容 10,469
关键字:

傅里叶变换与反变换的计算