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

编程问答 2022-05-15

问题描述

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

enter image description here

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,filtered_signals,ph,amp):
    equations = []
    for i in range(len(wavelength)):
        temp = eqn(filtered_signals[i],wavelength[i],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,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()

解决方法

这些是您的代码的问题：

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)
#                  ^^^

请注意，这两个更改纠正了绘制的图形，但并未纠正上述代码中的所有问题。

,

经过两天的挫折，我终于解决了这个问题。我仍然不知道为什么会这样，所以任何见解都会很棒。解决方法是使用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])

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

fftfftnumpyphasepythonsignal-processing