首页 > 学术百科

快速傅里叶变换及python代码实现

快速傅⾥叶变换及python代码实现

⼀、前⾔

我想认真写好快速傅⾥叶变换（Fast Fourier Transform，FFT），所以这篇⽂章会由浅到细，由窄到宽的讲解，但是傅⾥叶变换对于寻常⼈并不是很容易理解的，所以对于基础不牢的⼈我会通过前⾔普及⼀下相关知识。

我们复习⼀下三⾓函数的标准式：

y=A\cos (\omega x+\theta )+k

A代表振幅，函数周期是\frac{2\pi}{w}，频率是周期的倒数\frac{w}{2\pi}，\theta 是函数初相位，k在信号处理中称为直流分量。这个信号在频域就是⼀条竖线。

我们再来假设有⼀个⽐较复杂的时域函数y=f(t)，根据傅⾥叶的理论，任何⼀个周期函数可以被分解为⼀系列振幅A，频率\omega或初相位\theta 正弦函数的叠加

y = A_1sin(\omega_1t+\theta_1) + A_2sin(\omega_2t+\theta_2) + A_3sin(\omega_3t+\theta_3)

该信号在频域有三条竖线组成，⽽竖线图我们把它称为频谱图，⼤家可以通过下⾯的动画了解

如图可知，通过时域到频域的变换，我们得到了⼀个从侧⾯看的频谱，但是这个频谱并没有包含时域中全部的信息。因为频谱只代表每个正弦波对应频率的振幅是多少，⽽没有提到相位。基础的正弦波Asin(wt+\theta )中，振幅，频率，相位缺⼀不可，不同相位决定了波的位置，所以对于频域分析，仅仅有频谱（振幅谱）是不够的，我们还需要⼀个相位谱。

我依稀记得⾼中学正弦函数的是时候，\theta 的多少决定了正弦波向右移动多少。当然那个时候横坐标是相位⾓度，⽽时域信号的横坐标是时间，因此我们只需要将时间转换为相位⾓度就得到了初相位。相位差则是时间差在⼀个周期中所占的⽐例

\theta=2\pi \frac{t}{T}

所以傅⾥叶变换可以把⼀个⽐较复杂的函数转换为多个简单函数的叠加，将时域（即时间域）上的信号转变为频域（即频率域）上的信号，看问题的⾓度也从时间域转到了频率域，因此在时域中某些不好处理的地⽅，在频域就可以较为简单的处理，这就可以⼤量减少处理信号计算量。信号经过傅⾥叶变换后，可以得到频域的幅度谱以及相位谱，信号的幅度谱和相位谱是信号傅⾥叶变换后频谱的两个属性。

傅⾥叶⽤途

时域复杂的函数，在频域就是⼏条竖线

求解微分⽅程，傅⾥叶变换则可以让微分和积分在频域中变为乘法和除法

傅⾥叶变换相关函数

假设我们的输⼊信号的函数是

S=0.2+0.7*\cos (2\pi*50t+\frac{20}{180}\pi)+0.2*\cos (2\pi*100t+\frac{70}{180}\pi)

可以发现直流分量是0.2，以及两个余弦函数的叠加，余弦函数的幅值分别为0.7和0.2，频率分别为50和100，初相位分别为20度和70度。

freqs = np.fft.fftfreq(采样数量, 采样周期) 通过采样数与采样周期得到时域序列经过傅⾥叶变换后的频率序列

np.fft.fft(原序列) 原函数值的序列经过快速傅⾥叶变换得到⼀个复数数组，复数的模代表的是振幅，复数的辐⾓代表初相位

杀夫np.fft.ifft(复数序列) 复数数组经过逆向傅⾥叶变换得到合成的函数值数组

案例：针对合成波做快速傅⾥叶变换，得到分解波数组的频率、振幅、初相位数组，并绘制频域图像。

import matplotlib.pyplot as plt

import numpy as np

import numpy.fft as fft

Fs = 1000; # 采样频率

T = 1/Fs; # 采样周期

L = 1000; # 信号长度

t = [i*T for i in range(L)]

t = np.array(t)

S = 0.2+0.s(2*np.pi*50*t+20/180*np.pi) + 0.s(2*np.pi*100*t+70/180*np.pi) ;

complex_array = fft.fft(S)

nef

print(complex_array.shape) # (1000,)

print(complex_array.dtype) # complex128

print(complex_array[1]) # (-2.360174309695419e-14+2.3825789764340993e-13j)

>>>>>>###

plt.subplot(311)

plt.plot(1000*t[1:51], S[1:51], label='S') # y是1000个相加后的正弦序列

plt.xlabel("t（毫秒）")

plt.ylabel("S(t)幅值")

plt.title("叠加信号图")

plt.legend()

>>>>>>>

plt.subplot(312)

S_ifft = fft.ifft(complex_array)

# S_new是ifft变换后的序列

plt.plot(1000*t[1:51], S_ifft[1:51], label='S_ifft', color='orangered')

plt.xlabel("t（毫秒）")

plt.ylabel("S_ifft(t)幅值")

plt.title("ifft变换图")

plt.legend()

金惠敬

>>>>>>>

# 得到分解波的频率序列

freqs = fft.fftfreq(t.size, t[1] - t[0])

# 复数的模为信号的振幅（能量⼤⼩）

pows = np.abs(complex_array)

plt.subplot(313)

plt.title('FFT变换，频谱图')

plt.xlabel('Frequency 频率')

plt.ylabel('Power 功率')

plt.tick_params(labelsize=10)

plt.plot(freqs[freqs > 0], pows[freqs > 0], c='orangered', label='Frequency') plt.legend()

plt.tight_layout()

plt.show()

python代码实现

clear

clc

ADAMSGOLFclose all

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sampling period

L = 1000; % Length of signal

t = (0:L-1)*T; % Time vector

S = 0.2-0.7*cos(2*pi*50*t+20/180*pi) + 0.2*cos(2*pi*100*t+70/180*pi) ; plot(1000*t(1:50),S(1:50))

title('叠加信号图')

xlabel('t (milliseconds)')

ylabel('S(t)')

figure

Y = fft(S);

P2 = abs(Y/L);

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

plot(f,P1,'linewidth',2)

title('FFT变换')

xlabel('频率(Hz)')

ylabel('幅值')

figure

pred_X=ifft(Y);

plot(1000*t(1:50),pred_X(1:50),'r-')

MATLAB实现

补充⼀些复数知识（很重要）：

1、复数S的⼏种表⽰形式：

实部、虚部（直⾓坐标系）：a+bj （a是实部，b是虚部）

幅值、相位（指数系）：re^{j\theta } （r是幅值，\theta是相⾓，e^{j\theta }是相位）

极坐标表⽰法：r\angle \theta

指数系<-->指教坐标系：re^{j\theta }=r(cos\theta+jsin\theta)=rcos\theta+jrsin\theta

因此，我们可以通过以下⽅法得到：

实部：a=rcos\theta， real = np.real(S)

虚部：b=rsin\theta， imag= np.imag(S)

幅值：r=\sqrt{a^2+b^2}， magnitude = np.abs(S) 或 magnitude = np.sqrt(real**2+imag**2)

相⾓(以弧度为单位rad)：\theta=tan^{-1}(\frac{b}{a})或\theta=atan2(b,a)。 angle = np.angle(D(F, T))

相⾓(以⾓度为单位deg)：deg = rad*\frac{180}{\pi}，\text{rad2deg}(\text{atan2}(b,a))。 deg = rad * 180/np.pi

相位： phase = np.exp(1j * np.angle(S))

基于傅⾥叶变换的频域滤波

从某条曲线中除去⼀些特定的频率成份，这在⼯程上称为“滤波”。

含噪信号是⾼能信号与低能噪声叠加的信号，可以通过傅⾥叶变换的频域滤波实现降噪。

通过FFT使含噪信号转换为含噪频谱，去除低能噪声，留下⾼能频谱后再通过IFFT留下⾼能信号。

案例：基于傅⾥叶变换的频域滤波为⾳频⽂件去除噪声()。

1、读取⾳频⽂件，获取⾳频⽂件基本信息：采样个数，采样周期，与每个采样的声⾳信号值。绘制⾳频时域的：时间/位移图像import numpy as np

import numpy.fft as nf

import scipy.io.wavfile as wf

import matplotlib.pyplot as plt

# 读取⾳频⽂件

sample_rate, noised_sigs = wf.read('./da_data/noised.wav')

print(sample_rate) # sample_rate：采样率44100

print(noised_sigs.shape) # noised_sigs：存储⾳频中每个采样点的采样位移(220500,)

times = np.arange(noised_sigs.size) / sample_rate

plt.figure('Filter')

plt.subplot(221)

plt.title('Time Domain', fontsize=16)

plt.ylabel('Signal', fontsize=12)

plt.tick_params(labelsize=10)

plt.plot(times[:178], noised_sigs[:178], c='orangered', label='Noised')

plt.legend()

2、基于傅⾥叶变换，获取⾳频频域信息，绘制⾳频频域的：频率/能量图像

# 傅⾥叶变换后，绘制频域图像

grid servicefreqs = nf.fftfreq(times.size, times[1] - times[0])

complex_array = nf.fft(noised_sigs)

pows = np.abs(complex_array)

plt.subplot(222)

plt.title('Frequency Domain', fontsize=16)

plt.ylabel('Power', fontsize=12)

plt.tick_params(labelsize=10)

# 指数增长坐标画图

plt.semilogy(freqs[freqs > 0], pows[freqs > 0], c='limegreen', label='Noised')

plt.legend()

3、将低频噪声去除后绘制⾳频频域的：频率/能量图像

# 寻能量最⼤的频率值

fund_freq = freqs[pows.argmax()]

# where函数寻那些需要抹掉的复数的索引

noised_indices = np.where(freqs != fund_freq)

# 复制⼀个复数数组的副本，避免污染原始数据

filter_complex_array = py()

filter_complex_array[noised_indices] = 0

filter_pows = np.abs(filter_complex_array)

plt.subplot(224)

plt.xlabel('Frequency', fontsize=12)

plt.ylabel('Power', fontsize=12)

plt.tick_params(labelsize=10)

plt.plot(freqs[freqs >= 0], filter_pows[freqs >= 0], c='dodgerblue', label='Filter')

plt.legend()

4、基于逆向傅⾥叶变换，⽣成新的⾳频信号，绘制⾳频时域的：时间/位移图像filter_sigs = nf.ifft(filter_complex_array).real

plt.subplot(223)

plt.xlabel('Time', fontsize=12)

plt.ylabel('Signal', fontsize=12)

plt.tick_params(labelsize=10)

plt.plot(times[:178], filter_sigs[:178], c='hotpink', label='Filter')

律师事务所从事证券法律业务管理办法

plt.legend()

本文发布于:2024-09-23 20:14:41，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/436971.html

上一篇：连续时间傅里叶变换

下一篇：小波变换

标签：信号频域相位变换函数频率

留言与评论（共有 0 条评论）