References:
https://zhuanlan.zhihu.com/p/31584464
https://blog.csdn.net/weixin_39591031/article/details/110392352
https://blog.csdn.net/weixin_40146921/article/details/122006305
The NumPy module provides Fast Fourier Transform (FFT) functionality. FFT is an algorithm for efficiently calculating the discrete Fourier transform (DFT) and its inverse transform. It can quickly calculate the expression of time signals in the frequency domain. In NumPy, you can use the numpy.fft module to perform fast Fourier transform. This module provides a variety of functions for performing one-dimensional and multi-dimensional FFT operations. [For details, please refer to the numpy documentation]
np.fft.fft(a, n=None,axis=-1,norm=None)
This function computes the one-dimensional (n-element) discrete Fourier transform (DFT) using the Fast Fourier Transform algorithm
parameter:
- a: input array, which can contain complex numbers
- n: an integer, optional parameter. It represents the length of the final output of the selected axis. If n is smaller than len (the array of the axis of a), it will be clipped. If it is larger, fill it with 0s. n defaults to len (the array of the axis of a)
- axis: an integer, optional parameter. Indicates which axis of the input array a is selected. The last axis is selected by default.
- norm: normalization method, optional parameters. {“backward”, “ortho”, “forward”}, the default is “backward”, indicating that the backward direction of the transformation will be scaled
Output: n-element array with complex numbers of shape (1,n)
If x is a 1-dimensional array, then the FFT is equivalent to:
y[k] = np.sum(x * np.exp(-2j * np.pi * k * np.arange(n)/n))
-
y
k
=
∑
last pair
x
Sum each element of
x
?
e
?
2
π
j
?
k
?
[
0
,
1
,
.
.
.
,
n
?
1
]
n
y_k = \sum_{Finally sum each element of x} x \cdot e^{\frac{-2\pi j \cdot k \cdot [0,1,…,n -1]}{n}}
yk?=∑Finally, sum each element of x?x?en?2πj?k?[0,1,…,n?1]?
From the above equation we know the frequency
f
=
k
/
n
f=k/n
f=k/n
When k reaches
n
/
2
n/2
n/2 hours,
y
n
/
2
y_{n/2}
When yn/2?, it is already the Nyquist limit. At this time, the steering is performed, starting from k being the smallest negative number (
?
n
/
2
-n/2
?n/2) Start asking for help.
For example, if an x has 8 elements, then the calculation order of its corresponding k is as follows:
[0, 1, 2, 3, -4, -3, -2, -1]. You can also manually perform a translation fftshift from small to large. [-4, -3, -2, -1, 0, 1, 2, 3]
In this fft algorithm, a real number sequence x(n) is mapped to a complex number sequence X(n), and the output complex number sequence represents the frequency domain representation of the input real number sequence. Each complex number X(k) consists of a real part and an imaginary part. The imaginary part of the complex number represents some non-real components present in the input signal, such as some high-frequency noise or phase information.
Although the result of the FFT transform is a complex number, usually we only focus on the real part or amplitude information because this information is sufficient for many signal processing tasks. The imaginary part is usually ignored or used for some special applications such as phase recovery or phase calibration.
Usage example:
import numpy as np
import matplotlib
matplotlib.use('TkAgg')
matplotlib.rcParams['font.family'] = 'SimHei'
matplotlib.rcParams['font.size'] = 10
matplotlib.rcParams['axes.unicode_minus']=False
import matplotlib.pyplot as plt
# Draw the original signal
plt.figure(figsize=(8, 4))
plt.subplot(3, 1, 1)
plt.plot(np.linspace(0, 0.5, 200), np.cos(2 * np.pi * 5 * np.linspace(0, 0.5, 200)), label='original signal')
plt.xlabel('time')
plt.ylabel('magnitude')
plt.title('Original signal time domain representation')
plt.legend()
# Sample the original signal and perform FFT transformation
N = 20 # Sample 20 points
t = np.linspace(0, 0.5, N) # The time is 0-0.5s, there are N sample points, then the sampling frequency is (N-1) / 0.5 Hz
signal = np.cos(2 * np.pi * 5 * t) # The original signal is a cosine wave signal with frequency 5 and amplitude 1. Suppose we can collect N points of it
# Perform Fourier transform
fft_signal = np.fft.fft(signal)
# Calculate frequency axis
freq = np.fft.fftfreq(signal.shape[-1])
# freq = np.fft.fftfreq(signal.shape[-1], (t[-1] - t[0]) / (N - 1))
print(freq) # f = k / n, k = 0, 1, ... n/2, -n/2, -n/2 + 1, ...]
print(fft_signal)
# Draw frequency domain plot
plt.subplot(3, 1, 2)
plt.plot(np.fft.fftshift(freq), np.fft.fftshift(np.abs(fft_signal)), label='Frequency domain representation')
plt.xlabel('frequency')
plt.ylabel('magnitude')
plt.title('Frequency domain graph')
plt.legend()
signal_reconstructed = np.fft.ifft(fft_signal) # It is a complex number. When drawing, only the real part is taken.
print(signal_reconstructed)
plt.subplot(3, 1, 3)
plt.plot(t, np.real(signal_reconstructed), label='Reconstructed original signal')
plt.legend()
plt.xlabel('time')
plt.ylabel('magnitude')
plt.title('Reconstructed original signal time domain representation')
plt.tight_layout()
plt.show()
In the frequency domain diagram, when f is 0.15 and -0.15 (that is, when k is 3 and -3), the maximum amplitude is 8.51 (3.86-7.58j).
frequency formula
After knowing the sampling frequency Fs, the actual frequency corresponding to the xth (x starts from 0) complex value after fast Fourier transform (FFT) is
f(x) = x * (Fs / n)
For real signals, there is symmetry. Only half of the frequency 0 ~ N/2 is needed, and the other half is symmetrical with this half.
In other words, the actually calculated frequency is f = 0.15 * Fs= 0.15 * 38 = 5.7. Although this is not the actual frequency 5, it is very close. The main reason is that the sampling interval here is too small, resulting in low frequency resolution.
Let me briefly explain why: First of all, the sampling frequency is at least 5*2=10, which I am obviously satisfied with. But the time I took is 0-0.5s and N points, so that my sampling frequency is Fs = (N-1)/0.5=2N-1, and f(k ) = k * (Fs / N) = 2k-k/N, k is always rounded, so f(k)=5 will not occur at all, the closest one is f(3) =6-3/N, the key is to make the interval of 0.5 smaller.
The relationship between amplitude in the frequency domain and amplitude in the time domain is as follows:
Then the maximum amplitude after the above conversion to time domain should be 8.51 * (20 / 0.5) = 0.851
In other words, the frequency domain code can be modified as follows:
plt.plot(np.fft.fftshift(freq * ((N-1)/0.5) ), np.fft.fftshift(np.abs(fft_signal)/ (N/2) ), label=' Frequency domain representation') [Because there is no DC component here (there is no constant component in the expression), it is simply ignored]
Or use the parameter sampling time interval d when finding freq: freq = np.fft.fftfreq(signal.shape[-1], d= 1.0/ Fs)
So how to extract the main frequencies from the frequency domain after FFT transformation?
Here is a simple example of extracting the maximum frequency, using a periodic spectrum:
import numpy as np
import scipy
import scipy.io
from scipy.signal import butter
def select_max_freq(ori_signal, fs): # Select the largest freq in a small area in the frequency domain diagram
"""
using Fast Fourier transform (FFT).
"""
ori_signal = np.expand_dims(ori_signal, 0) # shape is converted from (n,) to (1, n) into a row vector
# It can be understood that the two-dimensional frequency domain map (frequency, amplitude) is compressed into a one-dimensional array [freq1_amplitude, freq2_amplitude, ...]
f, pxx = scipy.signal.periodogram(ori_signal, fs=fs, detrend=False)
# Select the frequency corresponding to the largest one from the power density
max_freq = np.take(f, np.argmax(pxx))
return max_freq
max_freq = select_max_freq(signal, Fs)
print(max_freq)
The final code is as follows:
import numpy as np
import matplotlib
import scipy
import scipy.io
from scipy.signal import butter
# matplotlib.use('TkAgg')
matplotlib.rcParams['font.family'] = 'SimHei'
matplotlib.rcParams['font.size'] = 10
matplotlib.rcParams['axes.unicode_minus'] = False
import matplotlib.pyplot as plt
# Draw the original signal
plt.figure(figsize=(8, 4))
plt.subplot(3, 1, 1)
plt.plot(np.linspace(0, 0.5, 200), np.cos(2 * np.pi * 5 * np.linspace(0, 0.5, 200)), label='original signal')
plt.xlabel('time')
plt.ylabel('magnitude')
plt.title('Original signal time domain representation')
plt.legend()
# Sample the original signal and perform FFT transformation
N = 20 # Sample 20 points
start_time = 0
stop_time = 0.5
t = np.linspace(start_time, stop_time, N) # The time is 0-0.5s, there are N sample points, the sampling frequency is (N-1) / 0.5 Hz
signal = np.cos(2 * np.pi * 5 * t) # The original signal is a cosine wave signal with frequency 5 and amplitude 1. Suppose we can collect N points of it
Fs = (N - 1) / (stop_time - start_time)
# Perform Fourier transform
fft_signal = np.fft.fft(signal)
# Calculate frequency axis
# freq = np.fft.fftfreq(signal.shape[-1])
freq = np.fft.fftfreq(signal.shape[-1], 1./Fs)
print(freq)
print(fft_signal)
# Draw frequency domain plot
plt.subplot(3, 1, 2)
plt.plot(np.fft.fftshift(freq), np.fft.fftshift(np.abs(fft_signal)), label='Frequency domain representation')
plt.xlabel('frequency')
plt.ylabel('magnitude')
plt.title('Frequency domain graph')
plt.legend()
signal_reconstructed = np.fft.ifft(fft_signal) # It is a complex number. When drawing, only the real part is taken.
# If you still want to use the imaginary part of the complex number, you can do this
#signal_reconstructed = np.piecewise(signal_reconstructed,
# [np.real(signal_reconstructed) < 0, np.real(signal_reconstructed) >= 0],
# [lambda elem: -np.abs(elem), lambda elem: np.abs(elem)]).astype(np.float64)
print(signal_reconstructed)
plt.subplot(3, 1, 3)
plt.plot(t, np.real(signal_reconstructed), label='Reconstructed original signal')
plt.legend()
plt.xlabel('time')
plt.ylabel('magnitude')
plt.title('Reconstructed original signal time domain representation')
plt.tight_layout()
plt.show()
def select_max_freq(ori_signal, fs): # Select the largest freq in a small area in the frequency domain diagram
"""
using Fast Fourier transform (FFT).
"""
ori_signal = np.expand_dims(ori_signal, 0) # shape is converted from (n,) to (1, n) into a row vector
# It can be understood that the two-dimensional frequency domain map (frequency, amplitude) is compressed into a one-dimensional array [freq1_amplitude, freq2_amplitude, ...]
f, pxx = scipy.signal.periodogram(ori_signal, fs=fs, detrend=False)
# Select the frequency corresponding to the largest one from the power density
max_freq = np.take(f, np.argmax(pxx))
return max_freq
max_freq = select_max_freq(signal, Fs)
print(max_freq) # 5.7