0、导入库

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import find_peaks

1、创建时域信号

创建时间序列

T = 0.01 # 采样间隔
fs = 100 # 采样频率
L = 1000 # 采样点数
tl = 0 # 起始时间
tr = tl + (L - 1) * T # 终止时间
t = np.linspace(tl, tr, L)

创建时域信号

以正弦信号为例

def generate_time_domain_sin_signal(t, A, f):pi = np.pix = A * np.sin(2 * pi * f * t)return x

A1 = 1.4
A2 = 0.6
f1 = 16
f2 = 17
x = generate_time_domain_sin_signal(t, A1, f1) + generate_time_domain_sin_signal(t, A2, f2)

画一下，看看效果

2、傅里叶变换

理论分析

• 原信号为正弦信号，假设频率为$f_0$的分量幅值为$A_0$，则傅里叶变换后，在f域中，在$f=\pm f_0$处的幅值应为$\frac{A_0}{2}$;相角在$f_0$处为$-90°$，在$-f_0$处为$+90°$
• 采样频率为$fs=100Hz$，相当于频域卷积梳状函数$fs{\delta }_{fs}(f)$，体现在幅值上，应乘$f_s$
• 原信号其实是无限信号，但是只采样了$L$个点，相当于在$t=L\ast T=10s$处截断了，相当于乘上了一个 门宽为10，门高为1的矩形窗，相当于频域卷积了一个峰值为10的sinc函数，幅值应乘$L\ast T$
• 综合以上，幅值应乘$\frac{L}{2}$
• 假如再在原信号的基础上加个别的窗，以hamming窗为例，幅值应乘hamming窗的积分，设积分为$I=\sum hamming[t]dt$，总幅值的变化倍数为$\frac{fs}{2}\ast I=\frac{\sum windowFunction[i]}{2}$

创建频率序列

在数字信号处理当中，任何信号的傅里叶变换都以奈奎斯特区间为周期，奈奎斯特区间为$[0,f_s]$，对应角频率为$w=[0,2pi]\ast fs$

批注：

奈奎斯特区间对应关系：

• $[0,1]$
• $[-0.5,0.5]$
• $[0,2\pi ]$
• $[-\pi ,\pi ]$
• $[0,fs]$
• $[-\frac{fs}{2},\frac{fs}{2}]$
• $[0,2\pi fs]$
• $[-\pi fs,\pi fs]$

w = np.linspace(0, 2 * pi, SAMPLE_N) * fs
f = w / 2 / pi

做DTFT

def DTFT(nT, xn, w):Xw = np.zeros(len(w), dtype=complex)for i, wi in enumerate(w):Xw[i] = np.sum(xn * np.exp(-1j * wi * nT))#  Xw = np.fft.fft(xn)return Xw
Xw = DTFT(t, x, w)

画一下，看看效果

幅值也满足我们的推导。

归一化

对结果乘$\frac{2}{L}$进行归一化

加窗

加个长度为$L$和hamming窗

def hamming_window(L):a0 = 0.53836return a0 - (1 - a0) * np.cos(2 * np.pi * np.arange(L) / (L - 1))

同时要准备好hanmming窗函数的积分，用于观察幅值变换

画个图，看看效果

时域

频域


非常nice

3、不同傅里叶变换的比较

1）CTFT

def CTFT(x, t, w):"""x[i] and t[i] is the i-th sample of the signal and time,for each w[i], this function will calculate the CTFT of x(t) at w[i].## Examples```pydef g(t):return np.where((t >= -4) & (t <= 4), 2, 0)t_values = np.linspace(-5, 9, SAMPLE_N)g_values = g(t_values)maxw = 10 * np.piw_values = np.linspace(-maxw, maxw, SAMPLE_N)Gw = CTFT(g_values, t_values, w_values)```"""Xw = np.zeros_like(w, dtype=complex)dt = t[1] - t[0]for i, wi in enumerate(w):# Two iterators here, x and tXw[i] = np.sum(x * np.exp(-1j * wi * t) * dt)return Xw

正弦谐波做完CTFT的结果：

相比于DTFT，CTFT在计算机当中的做法就是用DTFT实现CTFT的效果，但是多乘了个$dt$

2）DTFT

已做过

3）DFT

def DFT(ys):"""Calculate the Discrete Fourier Transform of the signal `ys`,which is an N-point sampling of the DTFT of the signal. ## Examples```pyy = np.array([-1, 2, 3, 0, -2, 1, 4, -3, 0, -2])dft_w_vec, dft_vec = DFT(y)debug_draw(dft_w_vec, np.abs(dft_vec))```"""n = len(ys)ns = np.arange(n)def omega_k(k):return 2 * np.pi * k / nw_vec = np.array([omega_k(k) for k in range(n)])dft_vec = np.array([sum(ys * np.exp(-1j * w * ns)) for w in w_vec])return w_vec, dft_vec

其中w_vec为数字角频率，范围为$[0,2pi]$，若要映射到范围$[0,fs]$,则需变换一下

w_dft, Xw_dft = DFT(x)
f_dft = w_dft / 2 / pi * fs

画一下，看看效果

同时我们还可以使用FFT实现

n_dft = len(x)
Xw_dft = np.fft.fft(x)
Xw_dft = np.fft.fftshift(Xw_dft)
f_dft = np.fft.fftfreq(n_dft, T)
f_dft = np.fft.fftshift(f_dft)

画一下，看看效果

4、帕斯瓦尔定理的验证

CTFT

print("帕斯瓦尔定理（CTFT）")
energy_time_domain = np.sum(np.abs(x) ** 2) * (t[1] - t[0])
energy_freq_domain = np.sum(np.abs(Xw) ** 2) * (f[1] - f[0])
print(energy_time_domain)
print(energy_freq_domain)
print("-------------")

DTFT

print("帕斯瓦尔定理（DTFT）")
energy_time_domain = np.sum(np.abs(x) ** 2)
energy_freq_domain = np.sum(np.abs(Xw) ** 2) / len(Xw)
print(energy_time_domain)
print(energy_freq_domain)
print("-------------")

DFT

print("帕斯瓦尔定理（DFT）")
energy_time_domain = np.sum(np.abs(x) ** 2)
energy_freq_domain = np.sum(np.abs(Xw_dft) ** 2) / n_dft
print(energy_time_domain)
print(energy_freq_domain)
print("-------------")

门函数的帕斯瓦尔定理

换一个信号，这次选择门函数，验证帕斯瓦尔定理是否正确

tl1 = 0
tr1 = 10
T = 0.01
fs = 100
L = 1000
t1 = np.linspace(tl1, tr1, L)
t1 = np.around(t1, 2)
x1 = np.where((t1 >= 2) & (t1 <= 8), 1, 0)w = np.linspace(0, 2 * pi, SAMPLE_N, endpoint=False) * fs

DTFT

Xw1 = DTFT(t1, x1, w)
print("帕斯瓦尔定理（DTFT）(门函数)")
energy_time_domain = np.sum(np.abs(x1) ** 2)
energy_freq_domain = np.sum(np.abs(Xw1) ** 2) / len(Xw1)
print(energy_time_domain)
print(energy_freq_domain)
print("-------------")

CTFT

Xw1 = CTFT(x1, t1, w)
print("帕斯瓦尔定理（CTFT）(门函数)")
energy_time_domain = np.sum(np.abs(x1) ** 2) * (t[1] - t[0])
energy_freq_domain = np.sum(np.abs(Xw1) ** 2) * (f[1] - f[0])
print(energy_time_domain)
print(energy_freq_domain)
print("-------------")

DFT

n1_dft = len(x1)
Xw1_dft = np.fft.fft(x1)
Xw1_dft = np.fft.fftshift(Xw1_dft)
f1_dft = np.fft.fftfreq(n1_dft, T)
f1_dft = np.fft.fftshift(f1_dft)
print("帕斯瓦尔定理（DFT）(门函数)")
energy_time_domain = np.sum(np.abs(x1) ** 2)
energy_freq_domain = np.sum(np.abs(Xw1_dft) ** 2) / n_dft
print(energy_time_domain)
print(energy_freq_domain)
print("-------------")

5、FFT

基本操作

对一个离散时间序列$x[n]$

SAMPLE_N = 5000
pi = np.pix = np.array([1, 2, -1, -3, 0, 2, 1, 1, 0, -1])
N = len(x)
n = np.arange(N)

做FFT，做完之后，通过fftshift将零频率移到数组的中间

Xw_dft = np.fft.fft(x)
Xw_dft = np.fft.fftshift(Xw_dft)
w_dft = np.fft.fftfreq(N, 1)
w_dft = np.fft.fftshift(w_dft)

比较DTFT与FFT的结果

zero-padding

在时间序列的最后进行zero_padding，可以增加点数且不改变频率信号的形状

x = np.append(x, np.zeros(20))
N = len(x)
n = np.arange(N)

再次画图

ifft

$np.fft.ifft$使用时需要注意：传进去的数组必须是[零频率项，正频率项，负频率项]
不能将fftshift过的数组传进去
画一下

跟原数组一样

zero-padding后也一样

6、FIR滤波器-频率采样法设计滤波器

创建时域信号

先准备一个时域信号，简单包含不同幅值，频率分别为$f1=5Hz,f2=5.5Hz$的两个正弦波叠加，采样频率为$fs=100Hz$，采$L=1000$个点

T = 0.01
fs = 100
A1 = 4
A2 = 2
f1 = 5
f2 = 5.5
L = 1000
t = np.linspace(0, (L - 1) * T, L)
t = np.around(t, 2)
x = generate_time_domain_sin_signal(t, A1, f1) + generate_time_domain_sin_signal(t, A2, f2)

画一画x[n]与X[f]

可以发现，频谱中有5Hz与5.5Hz两个分量

频域采样

设计一个低通滤波器，把5.5Hz滤掉。
在一个ideal的低通滤波器上采样，设置截止频率为5.25Hz，通带增益为2，阻带增益为0.1

def generate_lowpass_signal(fc, n_samples, passband, stopband):w = np.linspace(0, 2 * pi, n_samples)h = np.zeros(n_samples, dtype=complex)for i, wi in enumerate(w):if (wi <= fc or wi >= 2 * pi - fc):h[i] = passbandelse:h[i] = stopbandreturn h, w

orders = 1000
H, wh = generate_lowpass_signal(5.25 / (fs / 2) * pi, orders, 2, 0.1)

画一画

画出来也是挺ideal的！

反傅里叶变换

对采出来的滤波器数组进行$ifft$
要注意传进去的数组 遵循【零频率，正频率，负频率】顺序，我们是$[0,2\pi ]$，刚刚好！
做了$ifft$之后也要将零频率移到数组中间

h = np.fft.ifft(H)
h = np.fft.ifftshift(h)
n = np.arange(orders)

画一画，这里将阶数改为100了，因为1000太丑了

求滤波器的频率响应

先准备一个频率序列$w$，在$[-\pi ,\pi ]$中取几千个点

• 对于FIR滤波器，$h[n]=b[n],a[n]=[1,0,\mathrm{0...}]$
• 对于IIR滤波器，也不需要用频率采样哈哈哈

def freqz(bz, az, w):def cap_h_value(bz, az, z: complex):num = sum([bz[i] * z ** (-i) for i in range(len(bz))])den = sum([az[i] * z ** (-i) for i in range(len(az))])return num / dencap_hs = np.array([cap_h_value(bz, az, np.exp(1j * omega)) for omega in w])return cap_hs  # values of H(w)

w_frf = np.linspace(-pi, pi, SAMPLE_N)
H_frf = freqz(h, [1], w_frf)

画一画

发现还是频谱泄露了很多

缓解频谱泄漏：加个窗

对于一个滤波器序列，序列长度为$orders$，其零频率在序列中间
加个窗， 长度为$windowOrder$，通常加窗还有一个功能就是减少阶数，所以$windowOrder<orders$
巧的是，窗函数序列的零频率也在中间

所以有

def hamming_window(L):a0 = 0.53836return a0 - (1 - a0) * np.cos(2 * np.pi * np.arange(L) / (L - 1))

h = h[orders // 2 - window_order // 2 : orders // 2 + window_order // 2] * hamming_window(window_order)

这里以hamming窗为例了

加了窗之后再画一下频响，简直眉清目秀！

FIR滤波器的相位谱为线性相位，专门只取了200个点，看看相位：

滤波

def filter(bz, az, x, L):def expand(arr, n):return np.pad(arr, (0, n - len(arr)), 'constant')xs = expand(x, L)ys = np.zeros_like(xs)def y_at(n):return 0 if n < 0 else ys[n]def x_at(n):return 0 if n < 0 else xs[n]for n in range(L):ys[n] = sum([bz[i] * x_at(n - i) for i in range(len(bz))]) - \sum([az[i] * y_at(n - i) for i in range(1, len(az))])return xs, ys

x, y = filter(h, [1], x, L)

这里的阶数设置为500
滤波结果存在y里面

滤波结果分析

滤完之后，首先y的结果，需要把暂态的部分去掉，只看稳态部分

发现大概300~1000为稳态
把这一部分取出，再加窗缓解频谱泄漏

y = y[300 : 1000] * hamming_window(700)

画一画y的频谱

Yw = np.fft.fft(y)
Yw = np.fft.fftshift(Yw)
w = np.fft.fftfreq(700, T)
w = np.fft.fftshift(w)plot_a_signal(w, np.abs(Yw) * 2 / hamming_intergal(700))

发现5.5Hz的成分已经被消灭了，且幅值满足增益为2的设定。