我正在尝试理解傅里叶变换。我希望能够使用傅立叶变换提取时间序列数据集中的主导周期。
我几乎已经完成了,除了我无法获得函数的正确阶段,并且不知何故它们总是会移动。
在我的代码中我:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.fftpack import fft, ifft, fftfreq
import pandas as pd
# Generate a test signal
fs = 1000 # sampling frequency
time = np.linspace(0, 10, fs, False) # time vector
a_sine = np.sin(2 * np.pi * (time + 0.4) * 2)
b_sine = np.sin(2 * np.pi * (time + 0.2) * 0.6)
signal = a_sine + b_sine # test signal
###################################################################
N = len(signal)
f = fftfreq(N, 10/N) # i don't know why 10/N works but it gives the matching frequencies
signal_fft = fft(signal) # apply fourier transformation
# Get the magnitude and phase
magnitude = np.abs(signal_fft)
amplitude = magnitude / (N/2)
phase = np.angle(signal_fft)
dom, _ = find_peaks(amplitude, threshold=amplitude.max() * 0.3) # get the indices of dominant sine functions
dom = dom[:len(dom) // 2]
# get the parameters of sine functions
dominant_freqs = f[dom]
dominant_amplitude = amplitude[dom]
dominant_phases = phase[dom]
time = np.linspace(0, 10, N, endpoint=False)
# create sine functions
sine_functions = np.zeros_like(signal)
for freq, amp, phase in zip(dominant_freqs, dominant_amplitude, dominant_phases):
sine_functions += amp * np.sin(2 * np.pi * (time + phase) * freq)
# Created functions:
# 1 * np.sin(2 * np.pi * (time - 0.816) * 0.6)
# 1 * np.sin(2 * np.pi * (time - 2.827) * 2)
# as you can see extracted phases don't match the original ones: 0.4, 0.2
# Plot the original and reconstructed signals
plt.figure(figsize=(30, 15))
plt.subplot(2, 1, 1)
plt.plot(signal, label='Original Signal')
plt.plot(sine_functions, label='Reconstructed Signal')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(magnitude, label='Magnitude')
plt.show()
我现在将发布第一个更正,稍后继续。同样,采样频率的定义以及如何生成时间向量的问题可以通过以下方式解决:
# Generate a test signal
fs = 1000 # sampling frequency
TotalTime = 10 # total time
Ts = 1/fs # sampling time
time = np.linspace(0, TotalTime, TotalTime*fs, False) # time vector
a_sine = np.sin(2 * np.pi * (time + 0.4) * 2) # first sine
b_sine = np.sin(2 * np.pi * (time + 0.2) * 0.6) # second sine
signal = a_sine + b_sine # test signal
signalFFT = np.fft.fft(signal) # get fft
reconstructedSignal = np.fft.ifft(signalFFT) # reconstruct
plt.plot(signal) # plot signal
plt.plot(reconstructedSignal, "--") # reconstructed signal
plt.grid() # grid...
这是第一种重建方法
ifft
: