分解正弦波但相位有时似乎不匹配

Question

我最初在物理堆栈交换中发布了此内容，但他们要求也将其发布在这里.... 我正在尝试创建具有已知波长、幅度和相位的已知信号。然后我想把这个信号分解成它的所有频率，找到每个频率的幅度、相位和波长，然后根据这些新的波长、幅度和相位为每个频率创建方程。理论上，方程应该与单个信号相同。但是，它们不是。我几乎肯定这是一个阶段的问题，但我不知道如何解决它。我将发布确切的代码以在下面重现此内容。请帮忙，因为一旦我得到更复杂的信号，我的相位、波长和幅度会发生变化，因此它需要适用于这些的任何组合。

import numpy as np
from matplotlib import pyplot as plt
from scipy import fftpack

# create signal
time_vec = np.arange(1, 11, 1)
wavelength = 1/.1
phase = 0
amp = 10
created_signal = amp * np.sin((2 * np.pi / wavelength * time_vec) + phase)

# plot it
fig, axs = plt.subplots(2, 1, figsize=(10,6))
axs[0].plot(time_vec, created_signal, label='exact_data')

# get fft and freq array
sig_fft = fftpack.fft(created_signal)
sample_freq = fftpack.fftfreq(created_signal.size, d=1)

# do inverse fft and verify same curve as original signal. This is fine!
filtered_signal = fftpack.ifft(sig_fft)
filtered_signal += np.mean(created_signal)

# create individual signals for each frequency
filtered_signals = []
for i in range(len(sample_freq)):
    high_freq_fft = sig_fft.copy()
    high_freq_fft[np.abs(sample_freq) < np.nanmin(sample_freq[i])] = 0
    high_freq_fft[np.abs(sample_freq) > np.nanmax(sample_freq[i])] = 0
    filtered_sig = fftpack.ifft(high_freq_fft)
    filtered_sig += np.mean(created_signal)
    filtered_signals.append(filtered_sig)

# get phase, amplitude, and wavelength for each individual frequency
sig_size = len(created_signal)
wavelength = []
ph = []
amp = []
indices = []
for j in range(len(sample_freq)):
    wavelength.append(1 / sample_freq[j])
    indices.append(int(sig_size * sample_freq[j]))
for j in indices:
    phase = np.arctan2(sig_fft[j].imag, sig_fft[j].real)
    ph.append([phase])
    amp.append([np.sqrt((sig_fft[j].real * sig_fft[j].real) + (sig_fft[j].imag * sig_fft[j].imag)) / (sig_size / 2)])

# create an equation for each frequency based on each phase, amp, and wavelength found from above.
def eqn(filtered_si, wavelength, time_vec, phase, amp):
    return amp * np.sin((2 * np.pi / wavelength * time_vec) + phase)
def find_equations(filtered_signals_mean, high_freq_fft, wavelength, filtered_signals, time_vec, ph, amp):
    equations = []
    for i in range(len(wavelength)):
        temp = eqn(filtered_signals[i], wavelength[i], time_vec, ph[i], amp[i])
        equations.append(temp + filtered_signals_mean)
    return equations

filtered_signals_mean = np.abs(np.mean(filtered_signals))
equations = find_equations(filtered_signals_mean, sig_fft, wavelength, 
                                      filtered_signals, time_vec, ph, amp)

# at this point each equation, for each frequency should match identically each signal from each frequency,
# however, the phase seems wrong and they do not match!!??
axs[0].plot(time_vec, filtered_signal, '--', linewidth=3, label='filtered_sig_combined')
axs[1].plot(time_vec, filtered_signals[1], label='filtered_sig[-1]')
axs[1].plot(time_vec, equations[1], label='equations[-1]')
axs[0].legend()
axs[1].legend()
fig.tight_layout()
plt.show()  

Answer 1

这些是您的代码存在的问题：

filtered_signal = fftpack.ifft(sig_fft)
filtered_signal += np.mean(created_signal)

这仅适用于np.mean(created_signal) 近似为零。 IFFT 已经考虑了直流分量，零频率描述了信号的平均值。

filtered_signals = []
for i in range(len(sample_freq)):
   high_freq_fft = sig_fft.copy()
   high_freq_fft[np.abs(sample_freq) < np.nanmin(sample_freq[i])] = 0
   high_freq_fft[np.abs(sample_freq) > np.nanmax(sample_freq[i])] = 0
   filtered_sig = fftpack.ifft(high_freq_fft)
   filtered_sig += np.mean(created_signal)
   filtered_signals.append(filtered_sig)

在这里，在迭代的前半部分，遍历所有频率，同时考虑到负频率和正频率。例如，当i=1 时，您同时采用 -0.1 和 0.1 频率。您将 IFFT 应用于零信号的迭代的后半部分，根据定义，np.abs(sample_freq) 中没有一个小于零。

所以filtered_signals[1] 包含一个由 -0.1 和 0.1 频率分量构成的正弦波。这很好。否则它将是一个复值函数。

for j in range(len(sample_freq)):
   wavelength.append(1 / sample_freq[j])
   indices.append(int(sig_size * sample_freq[j]))

这里indices 数组的后半部分包含负值。不确定您对此有何计划，但它会导致后续代码从数组末尾开始索引。

for j in indices:
   phase = np.arctan2(sig_fft[j].imag, sig_fft[j].real)
   ph.append([phase])
   amp.append([np.sqrt((sig_fft[j].real * sig_fft[j].real) + (sig_fft[j].imag * sig_fft[j].imag)) / (sig_size / 2)])

这里，由于索引与前一个循环中的j 不同，phase[j] 并不总是对应于wavelength[j]，它们在大约一半的情况下引用来自不同频率分量的值。但我们不应该以任何方式评估这些案例。该代码假定一个实值输入，仅正频率的幅度和相位就足以重建信号。你应该跳过这里所有的负频率。

接下来，您使用收集的信息构建正弦波，但使用从 1 开始的 time_vec，而不是 FFT 假设的从 0 开始。因此，信号相对于预期值发生了偏移。此外，当phase==0 时，您应该创建一个偶数信号（即余弦，而不是正弦）。

因此，更改以下两行代码将创建正确的输出：

time_vec = np.arange(0, 10, 1)

和

def eqn(filtered_si, wavelength, time_vec, phase, amp):
   return amp * np.cos((2 * np.pi / wavelength * time_vec) + phase)
#                  ^^^

请注意，这两个更改更正了绘制的图形，但并没有更正上面讨论的代码中的所有问题。

Answer 2

经过 2 天的挫折，我终于解决了这个问题。我仍然不知道为什么会这样，所以任何洞察力都会很棒。解决方法是使用arctan2(Im, Re)产生的相位，并根据这个方程进行修改。

phase = np.arctan2(sig_fft[j].imag, sig_fft[j].real)    
formula = ((((wavelength[j]) / 2) - 2) * np.pi) / wavelength[j]
ph.append([phase + formula])

我不得不从数据中推导出这个方程，但我仍然不知道为什么会这样。请告诉我。终于！！