傅里叶变换理论基础，手撸源码

最新推荐文章于 2023-04-13 19:16:28 发布

weixin_38616018

最新推荐文章于 2023-04-13 19:16:28 发布

阅读量672

点赞数 2

分类专栏：编程图像文章标签：信号处理 python

本文链接：https://blog.csdn.net/weixin_38616018/article/details/107216682

版权

编程同时被 2 个专栏收录

8 篇文章 0 订阅

订阅专栏

图像

4 篇文章 0 订阅

订阅专栏

写在前边

我的人生中曾经两次接触过傅里叶，一次是我大学上高数的时候，那时候我赌考试不会考，于是根本没有看，事实证明我赌对了。另外一次是考研的时候，我记得很清楚是第7章无穷级数吧，里面有傅里叶级数，好像作用是把一个函数fx转换成无限个sinx和cosx的无穷级数吧，算是一个几百年没有考到过的知识点，我大致做过几道题，我还是赌他不会考，结果还是没有考。躲得过初一，躲不过十五，这次工作需要，我又拿起来课本。。。

数学基础

傅里叶变换所需要的数学概念有3个分别是正弦波，复数和点积，其中正弦波和复数可以延伸出复数正弦波而复数和正弦波可以延伸出复数点积这个概念，而这两个概念最终会让我们理解傅里叶变换

在这里插入图片描述

复数

首先说复数吧
我们知道对于一般的实数都是在一个一维的数轴上的
在这里插入图片描述我觉得复数就是就是数学家吃饱了撑的创造出来的一个集合，复数是没有任何意义的数字
就像x方什么时候等于-1一样，实际上就是无解，但是1维无解不代表2维没有解，就像我们永远无法突破四维空间一样，复数其实就是将1维的数轴变成了一个2维的平面，用一个一维不存在的i来表示另一个维度，虽然他在信号处理等领域很实用，扯远了，我相信大家对复数都比较清楚，毕竟大家都是经过9年义务教务的人哈哈。下图就是
在这里插入图片描述
复数可以让我们更有效的实现傅里叶变换，对于特定的频率，傅里叶变换会返回一个复数值，如下图

在这里插入图片描述
他到原点的距离表示振幅，他与x轴的夹角表示相位
代码,这边先引入各种后面要用的库

import numpy as np
import math
import matplotlib.pyplot as plt
import scipy.fftpack
import random
from mpl_toolkits.mplot3d import Axes3D
# import python library
z = complex(4,3)    # this way
print(z)
plt.plot(np.real(z),np.imag(z),'ro')
plt.axis('square')
plt.axis([-5, 5, -5, 5])
plt.grid(True)
plt.xlabel('Real axis'), plt.ylabel('Imaginary axis')
plt.show()
mag = np.abs(z)
print( '振幅',mag )
phs = np.angle(z)
print( '角度',phs )

欧拉公式

在这里插入图片描述
对于这个公式可以用这个图来表示

e的ix次方你可以想象成一个复数平面内长度为1方向为x的一个向量
对于这个欧拉公式其实我也想不太明白，我只是简单的验证了一下他

就当他是真理吧，毕竟涉及到了更高维（吃饱了撑）的复数了

x = np.linspace(-3,3,num=50)
plt.plot(x,np.exp(x),label='y=e^x')
plt.axis([min(x),max(x),0,np.exp(x[-1])])
plt.grid(True)
plt.legend()
plt.xlabel('x')
plt.show()

在这里插入图片描述

k = 2/np.pi
# Euler's notation
euler = np.exp(1j*k)
# plot dot
plt.plot(np.cos(k),np.sin(k),'ro')
# draw unit circle for reference
x = np.linspace(-np.pi,np.pi,num=100)
plt.plot(np.cos(x),np.sin(x))
# some plotting touch-ups
plt.axis('square')
plt.grid(True)
plt.xlabel('Real axis'), plt.ylabel('Imaginary axis')
plt.show()

在这里插入图片描述

正弦波和复数正弦波

正弦波相信大家都很熟悉了，毕竟9年义务教育，正弦波有频率和振幅还有相位，如果还是不懂那可以玩玩下面我写的代码看看改变那几个参数图像会变成啥样

# simulation parameters
srate = 500; # sampling rate in Hz
time  = np.arange(0.,2.,1./srate) # time in seconds
# sine wave param.eters
freq = 3;    # frequency in Hz
ampl = 2;    # amplitude in a.u.
phas = np.pi/3; # phase in radians
# generate the sine wave
sinewave = ampl * np.sin( 2*np.pi * freq * time + phas )
plt.plot(time,sinewave,'k')
plt.xlabel('Time (sec.)')
plt.ylabel('Amplitude (a.u.)')
plt.show()

需要注意的是频率和振幅还有相位是相互独立的，正弦和余弦是相互垂直的
那什么是复数正弦波呢？
很简单，我之前不是讲过欧拉公式吗
在这里插入图片描述
如果x用2πft+m带入
e的（2πft+m）i次方这个东西就是复数正弦波了

# general simulation parameters
srate = 500; # sampling rate in Hz
time  = np.arange(0.,2.,1./srate) # time in seconds
# sine wave parameters
freq = 5;    # frequency in Hz
ampl = 2;    # amplitude in a.u.
phas = np.pi/3; # phase in radians
# generate the sine wave
csw = ampl * np.exp( 1j* (2*np.pi * freq * time + phas) );
# plot the results
plt.plot(time,np.real(csw),label='real')
plt.plot(time,np.imag(csw),label='imag')
plt.xlabel('Time (sec.)'), plt.ylabel('Amplitude')
plt.title('Complex sine wave projections')
plt.legend()
plt.show()

在这里插入图片描述
这个就是复数正弦波了，这个代码执行后你可以看到实数部分和虚数部分随时间的变化，你可以试试把他变成3维然后画出来，哈哈我试了一下，是一个弹簧的形状。

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(time,np.real(csw),np.imag(csw))
ax.set_xlabel('Time (s)'), ax.set_ylabel('Real part'), ax.set_zlabel('Imag part')
ax.set_title('Complex sine wave in all its 3D glory')
plt.show()

在这里插入图片描述
那为什么要计算复数正弦波，如果我们用sin cos计算傅里叶变化的话，那我们的结果是取决于正弦波与信号之间相位的夹角的，如果我们稍微移动一下信号比如向后延迟5秒，尽管数据时一样的但是频谱图是完全不一样的

点积

点积我相信大部分人都知道，毕竟9年。。。
点积反应了一个相同维度下面两个向量的关系
如果点积大于0那么两个向量是锐角，小于0是钝角，等于0则是直角
那么点积的意义是什么呢
其实我们可以简单的理解成相关性
下面我们可以尝试一下一个信号和一个不同频率的实数正弦波的一个点积

# phase of signal
theta = 0*np.pi/4;
# simulation parameters
srate = 1000;
time  = np.arange(-1.,1.,1./srate)
# signal
sinew  = np.sin(2*np.pi*5*time + theta)
gauss  = np.exp( (-time**2) / .1);
signal = np.multiply(sinew,gauss)
# sine wave frequencies
sinefrex = np.arange(2.,10.,.5);
# plot signal
plt.plot(time,signal)
plt.xlabel('Time (sec.)'), plt.ylabel('Amplitude (a.u.)')
plt.title('Signal')
plt.show()
# initialize dot products vector
dps = np.zeros(len(sinefrex));
# loop over sine waves
for fi in range(1,len(dps)):  
    # create sine wave
    sinew = np.sin( 2*np.pi*sinefrex[fi]*time)   
    # compute dot product
    dps[fi] = np.dot( sinew,signal ) / len(time)
# and plot
plt.stem(sinefrex,dps)
plt.xlabel('Sine wave frequency (Hz)'), plt.ylabel('Dot product')
plt.title('Dot products with sine waves')
plt.show()

在这里插入图片描述
如果我们去改信号的相位也就是变量theta的值我们会发现输出的图片完全不一样了，所以这是有问题的，我们如果只用一个实数正弦波去做点积的话，他是非常取决于信号的相位的，相位发生了改变，点积会发生很大的改变，那就会导致我们的结果并不能表示信号的特征，那怎么解决这个问题呢，我们引入了复数点积

复数点积

复数点积是什么？
复数点积其实就是对信号的实数部分和虚数部分分别做点积

在这里插入图片描述
他到原点的距离表示振幅，同时也表示点积

# phase of signal
theta = 0*np.pi/4;
# simulation parameters
srate = 1000;
time  = np.arange(-1.,1.,1./srate)
# signal
sinew  = np.sin(2*np.pi*5*time + theta)
gauss  = np.exp( (-time**2) / .1);
signal = np.multiply(sinew,gauss)
# sine wave frequencies
sinefrex = np.arange(2.,10.,.5);
# plot signal
plt.plot(time,signal)
plt.xlabel('Time (sec.)'), plt.ylabel('Amplitude (a.u.)')
plt.title('Signal')
plt.show()
# initialize dot products vector
dps = np.zeros(len(sinefrex));
# loop over sine waves
for fi in range(1,len(dps)):  
    # create sine wave
    sinew = np.exp( 1j*2*np.pi*sinefrex[fi]*time )   
    # compute dot product
    dps[fi] = np.abs( np.dot( sinew,signal ) / len(time) )
# and plot
plt.stem(sinefrex,dps)
plt.xlabel('Sine wave frequency (Hz)'), plt.ylabel('Dot product')
plt.title('Dot products with sine waves')
plt.show()

在这里插入图片描述
这个代码其实和上面几乎一样，只是我把正弦波改成了复数正弦波，这时我们可以试着像上面一样去改变他的代码，我们发现无论我们怎么改信号的相位也就是变量theta，点积的结果没有任何变化
所以为什么会这样呢
这是因为复数正弦波有两个波，实数和虚数，当相位发生变化时，他的点永远是复平面上的一个点，相位表示这个向量与实数轴的一个夹角，但是他永远会在这个复平面上，现在我们算的是这个点到原点的距离，他的幅值是不会发生变化的

# create complex sine wave
csw = np.exp( 1j*2*np.pi*5*time )
rsw = np.sin(    2*np.pi*5*time )
# specify range of phase offsets for signal
phases = np.linspace(0,7*np.pi/2,num=100)
for phi in range(0,len(phases)):   
    # create signal
    sinew  = np.sin(2*np.pi*5*time + phases[phi])
    gauss  = np.exp( (-time**2) / .1)
    signal = np.multiply(sinew,gauss)
    # compute complex dot product
    cdp = np.sum( np.multiply(signal,csw) ) / len(time)
    # compute real-valued dot product
    rdp = sum( np.multiply(signal,rsw) ) / len(time)
    # plot signal and real part of sine wave
    pl.cla() # wipe the figure
    plt.subplot2grid((2,2), (0, 0), colspan=2)
    plt.plot(time,signal)
    plt.plot(time,rsw)
    plt.title('Signal and sine wave over time')
    # plot complex dot product
    plt.subplot2grid((2,2), (1, 0))
    plt.plot(np.real(cdp),np.imag(cdp),'ro')
    plt.xlabel('Real')
    plt.ylabel('Imaginary')
    plt.axis('square')
    plt.grid(True)
    plt.axis([-.2,.2,-.2,.2])
    plt.plot([0,np.real(cdp)],[0,np.imag(cdp)],'r')
    # draw normal dot product
    plt.subplot2grid((2,2), (1, 1))
    plt.plot(rdp,0,'ro')
    plt.xlabel('Real')
    plt.axis('square')
    plt.grid(True)
    plt.axis([-.2,.2,-.2,.2])
    # show plot    
    display.clear_output(wait=True)
    display.display(pl.gcf())
    ttime.sleep(.01)

运行这行代码我们可以发现信号的移动仅仅只改变右下角的实数部分图，对于左小角的图来说，它对应的虚数点是一个复平面上的一个圆，相位对应的是他的角度，而是到原点距离也就是他的点积是不会发生改变的。到这里我们已经掌握所有的数学基础了。

开始手撸傅里叶变换

傅里叶变换到底是怎么工作的呢
本来想直接翻译opencv里面的dft源码的，但是由于opencv的设计问题，我看了两天放弃了，实在看不懂啊哈哈。

我这里直接给出算法的伪代码吧

对于每个时间点的循环开始
创建一个与信号时间长度一样的复数正弦波，其中频率由信号点的个数减1决定
计算复数正弦波和信号之间的点积
结束
信号的振幅是傅里叶系数的幅值
信号的相位是傅里叶系数的角度

# create the signal
srate  = 100 # hz
time   = np.arange(0.,2.,1/srate) # time vector in seconds
pnts   = len(time) # number of time points
signal = 2 * np.sin( 2*np.pi*2*time ) + 4 * np.sin( 2*np.pi*4*time )
# prepare the Fourier transform
fourTime = np.array(range(0,pnts))/pnts
fCoefs   = np.zeros((len(signal)),dtype=complex)
for fi in range(0,pnts):  
    # create complex sine wave
    csw = np.exp( -1j*2*np.pi*fi*fourTime )  
    # compute dot product between sine wave and signal
    # these are called the Fourier coefficients
    fCoefs[fi] =  np.multiply(signal,csw).sum()  / pnts
# extract amplitudes
ampls = 2*np.abs(fCoefs)
# compute frequencies vector
hz = np.linspace(0,srate/2,num=math.floor(pnts/2.)+1)
plt.stem(hz,ampls[range(0,len(hz))])
plt.xlabel('Frequency (Hz)'), plt.ylabel('Amplitude (a.u.)')
plt.xlim(0,10)
plt.show()

这里我用numpy的fft函数做验证，可以发现上下两段代码输出的图是一样的

fCoefs = np.fft.fft(signal)
# extract amplitudes
ampls = 2*np.abs(fCoefs)/pnts
# compute frequencies vector
hz = np.linspace(0,srate/2,num=math.floor(pnts/2.)+1)
plt.stem(hz,ampls[range(0,len(hz))])
plt.xlabel('Frequency (Hz)'), plt.ylabel('Amplitude (a.u.)')
plt.xlim(0,10)
plt.show()

在这里插入图片描述
这边其实还有很多细节，比如为什么np.abs(fCoefs)要乘以2呢？
这是因为振幅不仅去了正频率也去了负频率，而正频率和负频率的关系其实是个镜像
还有很多细节涉及到“奈奎斯特频率”和“正规化因子”这里由于篇幅有限，就不展开了

快速傅里叶变换

上面我实现了离散型傅里叶变换DFT，但是，这种基于循环的傅里叶变换，面对大量的数据是力不从心的。我们实际使用，往往用的是快速傅里叶变换，是以稀疏矩阵为主要的数据结构，时间空间都更加合理。但是其实这两者的本质是一样的

对于二维图像

对于二维图像做傅里叶变化其实是对每一行做一次傅里叶之后再对每一列做一次傅里叶变换，然后一般低频放中间高频放两边，然后在log这样处理一下。其中低频信号表示大致的一个轮廓，而高频信号表示细节。

weixin_38616018

关注

2
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
傅里叶变换理论基础，手撸源码

写在前边我的人生中曾经两次接触过傅里叶，一次是我大学上高数的时候，那时候我赌考试不会考，于是根本没有看，事实证明我赌对了。另外一次是考研的时候，我记得很清楚是第7章无穷级数吧，里面有傅里叶级数，好像作用是把一个函数fx转换成无限个sinx和cosx的无穷级数吧，算是一个几百年没有考到过的知识点，我大致做过几道题，我还是赌他不会考，结果还是没有考。躲得过初一，躲不过十五，这次工作需要，我又拿起来课本。。。...
复制链接

扫一扫