I am trying to create a known signal with a known wavelength, amplitude, and phase. I then want to break this signal apart into all of its frequencies, find amplitudes, phases, and wavelengths for each frequency, then create equations for each frequency based on these new wavelengths, amplitudes, and phases. In theory, the equations should be identical to the individual signals. However, they are not. I am almost positive it is an issue with phase but I cannot figure out how to resolve it. I will post the exact code to reproduce this below. Please help as my phase, wavelength, and amplitudes will vary once I get more complicated signals so it need to work for any combination of these.

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)

# 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.angle(sig_fft[j])
    phase = np.arctan2(sig_fft[j].imag, sig_fft[j].real)
    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]')
  • 2
    $\begingroup$ arctan usually returns a number in the range $-{\pi}/2$ to ${\pi}/2$. Do you have an arctan2 available? arctan2(Re,Im) should fix the range ambiguity. (I haven't gone through all the details you posted, but this is a common problem in the ballpark of what you're talking about. If this isn't it, please mention it.) $\endgroup$ – tom10 Apr 3 at 4:39
  • 1
    $\begingroup$ To complement @tom10's comment, see my answer here that discusses why one should use the $\texttt{atan2}$ function in some detail. $\endgroup$ – Philip Apr 3 at 5:43
  • $\begingroup$ atan2 gives identical result. I also have a built in function that does this calculation in python. it is numpy.angle. All three of these give the same result. $\endgroup$ – Justin Oberle Apr 3 at 15:30
  • $\begingroup$ Is there a way to move this to the signal processing stack exchange, or should I just make a new one there? $\endgroup$ – Justin Oberle Apr 3 at 15:33
  • $\begingroup$ Ok. I will post in there. In the mean time, I broke down my code significantly to only the required pieces for duplication. I am going to edit the above post and show the code here. I will also add this to stackoverflow. If I don't hear anything over there, I will let you know. Thanks $\endgroup$ – Justin Oberle Apr 4 at 17:37

I solved this finally after 2 days of frustration. I still have no idea why this is the way it is so any insight would be great. The solution is to use the phase produced by arctan2(Im, Re) and modify it according to this equation.

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

I had to derive this equation from data but I still do not know why this works. Please let me know. Finally!!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.