AI专业者开发社区 1 号成员
  •  0 回帖  •  11 浏览  •  4 个月前
   

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

目录

一、前言

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

　　我们复习一下三角函数的标准式：

y=Acos⁡(ωx+θ)+k

　　A 代表振幅，函数周期是 2πw，频率是周期的倒数 w2π，θ是函数初相位，k 在信号处理中称为直流分量。这个信号在频域就是一条竖线。

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

y=A1sin(ω1t+θ1)+A2sin(ω2t+θ2)+A3sin(ω3t+θ3)

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

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

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

θ=2πtT

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

傅里叶用途

  • 时域复杂的函数，在频域就是几条竖线
  • 求解微分方程，傅里叶变换则可以让微分和积分在频域中变为乘法和除法

傅里叶变换相关函数

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

S=0.2+0.7∗cos⁡(2π∗50t+20180π)+0.2∗cos⁡(2π∗100t+70180π)

可以发现直流分量是 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

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示符号

Fs = 1000;            # 采样频率
T = 1/Fs;             # 采样周期
L = 1000;             # 信号长度
t = [i*T for i in range(L)]
t = np.array(t)


S = 0.2+0.7np.cos(2np.pi50t+20/180np.pi) + 0.2np.cos(2np.pi100t+70/180np.pi) ;


complex_array = fft.fft(S)
print(complex_array.shape)  # (1000,) 
print(complex_array.dtype)  # complex128 
print(complex_array[1])  # (-2.360174309695419e-14+2.3825789764340993e-13j)


#################################
plt.subplot(311)
plt.grid(linestyle=':')
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.grid(linestyle=':')
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.grid(linestyle=':')
plt.plot(freqs[freqs > 0], pows[freqs > 0], c='orangered', label='Frequency')
plt.legend()
plt.tight_layout()
plt.show()

python 代码实现

clear
clc
close 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) = 2P1(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 是虚部）
  • 幅值、相位（指数系）：rejθ　　（r 是幅值，θ是相角，ejθ是相位）
  • 极坐标表示法：r∠θ
  • 指数系 <−−> 指教坐标系：rejθ=r(cosθ+jsinθ)=rcosθ+jrsinθ

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

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

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

幅值：r=a2+b2， magnitude = np.abs(S)  或  magnitude = np.sqrt(real**2+imag**2) 

相角(以弧度为单位 rad)：θ=tan−1(ba) 或 θ=atan2(b,a)。 angle = np.angle(D(F, T)) 

相角(以角度为单位 deg)：deg=rad∗180π，rad2deg(atan2(b,a))。 deg = rad * 180/np.pi  

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

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

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

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

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

案例：基于傅里叶变换的频域滤波为音频文件去除噪声 (noiseed.wav 数据集地址)。

　　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.grid(linestyle=':')
plt.plot(times[:178], noised_sigs[:178], c='orangered', label='Noised')
plt.legend()

 

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

# 傅里叶变换后，绘制频域图像
freqs = 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.grid(linestyle=':')
# 指数增长坐标画图
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 = complex_array.copy()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.grid(linestyle=':')
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.grid(linestyle=':')
plt.plot(times[:178], filter_sigs[:178], c='hotpink', label='Filter')plt.legend()

　　5、重新生成音频文件

# 生成音频文件
wf.write('./da_data/filter.wav', sample_rate, filter_sigs)plt.show()

 

 

离散傅里叶变换 (DFT)

　　离散傅里叶变换（DFT）对有限长时域离散信号的频谱进行等间隔采样，频域函数被离散化了，便于信号的计算机处理。DFT 的运算量太大，FFT 是离散傅里叶变换的快速算法。

　　说白了 FFT 和 DFT 它俩就是一个东东，只不过复杂度不同，

　　有时候我们能够看到 N 点傅里叶变换，我个人理解是这个 N 点是信号前面 N 个连续的数值，即 N 点 FFT 意思就是截取前面 N 个信号进行 FFT，这样就要求我们的前 N 个采样点必须包含当前信号的一个周期，不然提取的余弦波参数与正确的叠加波的参数相差很大。

　　如果在 N 点 FFT 的时候，如果这 N 个采样点不包含一个周期呢？或者说我们的信号压根不是一个周期函数咋办？或者有一段是噪音数据呢？如果用 FFT 计算，就会对整体结果影响很大，然后就有人想通过局部来逼近整体，跟微积分的思想很像，将信号分成一小段一小段，然后对每一小段做 FFT，就跟分段函数似的，无数个分段函数能逼近任意的曲线 ((⊙o⊙)…应该没错吧 )，这样每一段都不会互相影响到了。

二、短时傅里叶变换 (STFT)

　　在短时傅里叶变换过程中，窗的长度决定频谱图的时间分辨率和频率分辨率，窗长越长，截取的信号越长，信号越长，傅里叶变换后频率分辨率越高，时间分辨率越差；相反，窗长越短，截取的信号就越短，频率分辨率越差，时间分辨率越好，也就是说短时傅里叶变换中，时间分辨率和频率分辨率之间不能兼得，应该根据具体需求进行取舍。

计算短时傅里叶变换，需要指定的有：

  • 每个窗口的长度：nsc
  • 每相邻两个窗口的重叠率：nov
  • 每个窗口的 FFT 采样点数：nff

可以计算的有：

  • 海明窗：w=hamming(nsc, 'periodic')
  • 信号被分成了多少片
  • 短时傅里叶变换：

python 库librosa实现

librosa.stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann', center=True, pad_mode='reflect')

短时傅立叶变换（STFT），返回一个复数矩阵使得 D(f,t)

参数：

  • y：音频时间序列
  • n_fft：FFT 窗口大小，n_fft=hop_length+overlapping
  • hop_length：帧移，如果未指定，则默认 win_length / 4。
  • win_length：每一帧音频都由 window（）加窗。窗长 win_length，然后用零填充以匹配 N_FFT。默认win_length=n_fft。
  • window：字符串，元组，数字，函数 shape =（n_fft, )
    • 窗口（字符串，元组或数字）；
    • 窗函数，例如scipy.signal.hanning
    • 长度为 n_fft 的向量或数组
  • center：bool
    • 如果为 True，则填充信号 y，以使帧 D [:, t] 以 y [t * hop_length] 为中心。
    • 如果为 False，则 D [:, t] 从 y [t * hop_length] 开始
  • dtype：D 的复数值类型。默认值为 64-bit complex 复数
  • pad_mode：如果 center = True，则在信号的边缘使用填充模式。默认情况下，STFT 使用 reflection padding。

返回：

  • STFT 矩阵 D，shape =（1 + nfft2，t）
librosa.magphase(D, power=1)

将复值频谱 D 分离成其幅度（S）和相位（P）

参数：

  • D：复值谱图，np.ndarray [shape =（d，t），dtype = complex]
  • power：幅度谱图的指数，例如，1 代表能量，2 代表功率，等等。

返回：

  • D_mag ：D 的幅值，np.ndarray [shape =（d，t），dtype = real]
  • D_phase ：D 的相位，np.ndarray [shape =（d，t），dtype = complex]，exp（1.j * phi）其中phi是D的相位
y, _ = librosa.load("p225_002.wav", sr=16000)

D = librosa.stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann', center=True, pad_mode='reflect')
print(D.shape)    # (1025, 127)


# 将复值频谱 D 分离成其幅度 (S) 和相位 (P) 的部件
magnitude, phase = librosa.magphase(D, power=1)
# magnitude        # 赋值 [shape =（d，t），dtype = real] (1025, 127)
# phase.shape    # 相位 [shape =（d，t），dtype = complex] (1025, 127) 

angle = np.angle(phase)     # 获取相角（以弧度为单位）

 

tensorflow 实现

一句话实现分帧、加窗、STFT 变换

# [batch_size, signal_length].  batch_size and signal_length 可能会不知道
signals = tf.placeholder(tf.float32, [None, None])

# stfts 短时傅里叶变换：就是对信号中每一帧信号进行傅里叶变换 
# shape is [batch_size, ?, fft_unique_bins]
# 其中 fft_unique_bins = fft_length // 2 + 1 = 513.
stfts = tf.contrib.signal.stft(signals, frame_length=1024, frame_step=512,
fft_length=1024)

wlen：窗长

frame_length 是信号中帧的长度

frame_step 是帧移

fft_length：做 fft 变换的长度，或一种说话：fft 变换所取得 N 点数，在有些地方也表示为 NFFT。

注意：FFT 的长度必须大于或者等于 win 的长度或者帧长。以获得更高的频域分辨率

FFT 后的分辨率（频率的间隔）为 fs/NFFT。当 NFFT>wlen 时就是在数据补零后做 FFT，做的 FFT 得到的频谱等于以 wlen 长数据 FFT 的频谱中内插。

numpy 库实现

上面的一行代码相当于下面一大段代码

def wav_to_frame(wave_data, win_len, win_shift):
    """
    进行分帧操作
    :param wave_data: 原始的数据
    :param win_len: 滑动窗长
    :param win_shift: 滑动间隔
    :return: 分帧之后的结果，输出一个帧矩阵
    """
    num_frames = (len(wave_data) - win_len) // win_shift + 1
    results = []
    for i in range(num_frames):
        results.append(wave_data[i*win_shift:i*win_shift + win_len])
    return np.array(results)

def spectrum_power(frames, NFFT):
"""
计算每一帧傅立叶变换以后的功率谱
参数说明：
frames:audio2frame 函数计算出来的帧矩阵
NFFT:FFT 的大小
"""
# 功率谱等于每一点的幅度平方 /NFFT
return 1.0/NFFT * np.square(spectrum_magnitude(frames, NFFT))


def spectrum_magnitude(frames, NFFT):
"""计算每一帧经过 FFT 变幻以后的频谱的幅度，若 frames 的大小为 NL, 则返回矩阵的大小为 NNFFT
参数：
frames: 即 audio2frame 函数中的返回值矩阵，帧矩阵
NFFT:FFT 变换的数组大小, 如果帧长度小于 NFFT，则帧的其余部分用 0 填充铺满
"""
complex_spectrum = np.fft.rfft(frames, NFFT)    # 对 frames 进行 FFT 变换
# 返回频谱的幅度值
return np.absolute(complex_spectrum)

View Code

 

三、frequency bin

在读 paper 的时候总是会遇到 frequency bin （频率窗口）这个词，

　　frequency bin 是指：raw data 经过 FFT 后得到的频谱图（频域率）中，频率轴的频率间隔或分辨率，通常取决采样率和采样点。

采样率采样点数采样率采样点数 frequencybin= 采样率采样点数 =fsampleNrecode

　　Nrecode  是信号在时域的采样点数，频谱中的频率点或线的数量为 Nrecode2  （奈奎斯特采样定理）

　　频谱的第一个频点始终为直流（频率 =0），最后一个频点为 fsample2−fsampleNrecode  。频点采用相等的间隔，这间隔通常用frequency bin（频率窗口）或 FFT bin 表示。

例子 1：我们可以作用 82MHz 的采样频率，取得 8200 个数据记录，frequency bin=820000008200=10000=10kHz。

例子 2：frequency bin 是频率域中采样点之间的间隔。例如，如果采样率为 100 赫兹，FFT 为 100 个点，frequency bin=1，则在 [0 100）赫兹之间有 100 个点。因此，您将整个 100 赫兹范围划分为 100 个间隔，如 0-1 赫兹、1-2 赫兹等。每一个如此小的间隔，比如 0-1Hz，都是一个 frequency bin（频率箱）。

参考

知乎 Heinrich 博主文章——傅里叶分析之掐死教程

【音频处理】离散傅里叶变换

【音频处理】短时傅里叶变换

作者：凌逆战
欢迎任何形式的转载，但请务必注明出处。
限于本人水平，如果文章和代码有表述不当之处，还请不吝赐教。
本文章不做任何商业用途，仅作为自学所用，文章后面会有参考链接，我可能会复制原作者的话，如果介意，我会修改或者删除。