1

Hello, I am manually compute the inverse FFT using the NumPy FFT components by the below codes, to compare the rebuild signal with the original signal (sin function). As shown in the results, the error (original - rebuilt) is increasing along the x axis. Does anyone know the reasons and how to correct it? If I used the build-in function ifft to rebuild the signal, the error doesn't have this increase tread along x axis. 

To show the error clearly, here I used sin signal as a simple example.my real signal is much more complicated. Also, I like to emphases I have to use the cos or exp analytical equation to rebuild the signal because I need this equation for further calculation, that is to say I can't use build in ifft function.

Although the absolute error value is pretty small, the increasing tread will make the results unusable for my application.  

Greatly Appreciate your help

import numpy as np
from scipy import *
from matplotlib import pyplot as plt
tempdata = np.loadtxt("sinx_80.dat")
x_positions = tempdata[:, 0]
delta_x=x_positions[1] - x_positions[0]
roughness_data = tempdata[:, 1]
fft_1d = np.fft.fft(roughness_data)
num_x = len(fft_1d)

freq_x = np.fft.fftfreq(num_x, d=delta_x)
reconstructed_surface=np.zeros(num_x)
x_grid = (x_positions - x_positions[0])

for n in range(num_x):
    amplitude = np.abs(fft_1d[n]) / (num_x)
    phase = np.angle(fft_1d[n])
    reconstructed_surface += amplitude * np.cos(2 * np.pi * (freq_x[n] * x_grid) + 
phase)

difference = roughness_data - reconstructed_surface
max_error = np.max(np.abs(difference))
rms_error = np.sqrt(np.mean(difference**2))

print(f"Max Error: {max_error}")
print(f"RMS Error: {rms_error}")
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(x_positions, roughness_data, label="Original Data", color="blue")
plt.title("Original Data")
plt.legend()
plt.grid()
plt.subplot(3, 1, 2)
plt.plot(x_positions, reconstructed_surface, label="Reconstructed Data", 
color="orange")
plt.title("Reconstructed Data")
plt.legend()
plt.grid()
plt.subplot(3, 1, 3)
plt.plot(x_positions, difference, label="Difference (Original - Reconstructed)", 
color="red")
plt.title("Difference")
plt.legend()
plt.grid()
plt.show()

Hi liwang,

The reconstruction formula is wrong.  Since a cosine function is being used as the basis, the sum should be performed only over the DFT coefficients corresponding to the positive frequencies (i.e., from 0 Hz up to the Nyquist frequency). In addition, the correct scaling factor for the amplitude is num_x/2. The resulting error, although not zero, should not increase with x.

Saludos

Saludos, Thank you so much for your help. Following your suggestion, I adjusted my codes, but unfortunately, the results are the same (first group plot).   For comparison, if using the buidin numpy function IFFT to rebuild the signal, this increasing tread disappeared (second group plot). unfortunately, I have to use analytical equation(sin, cos or exp) for further calculation.  Do you have insight for the reasons of increasing tread coming with the analytical equation? 

fft_1d = np.fft.fft(roughness_data)
num_x = len(fft_1d)
freq_x = np.fft.fftfreq(num_x, d=delta_x)
nyquist_freq = 1 / (2 * delta_x)
x_grid = x_positions - x_positions[0]
reconstructed_surface = np.zeros(num_x)
# Manually rebuild by only POSITIVE frequencies within the Nyquist limit
for n in range(num_x):
    if freq_x[n] > 0 and freq_x[n] <= nyquist_freq:  # Only positive within Nyquist limit
        amplitude = 2 * np.abs(fft_1d[n]) / num_x  # Multiply by 2 to compensate for ignoring negative frequencies
        phase = np.angle(fft_1d[n])
        reconstructed_surface += amplitude * np.cos(2 * np.pi * freq_x[n] * x_grid +phase)

Dear liwang, as pointed out before, the result will not generate an error strictly equal to zero. In addition, from your last error plot, it is not increasing with x, it remains bounded to +/-5e-20. Yo should see that this bounded behavior persist for longer x, in comparison to your initial code.

Hi ing_jpu, thanks for your reply. 

unfortunately the last error plot(which does not increase with x) is from buidin IFFT inverse function, and it doesn't fit my calculation purpose. I need analytical equation(sin, cos or exp) of inverse FFT function for further calculation.

The corrected implementation should not generate an increasing error. If it does, you should debug your code.

Suerte

I did not read your code but running a sin/cos function for very large argument can be increasingly inaccurate.

Knutlnge, appreciate your help. I also tried exp instead of sin/cos with little improvement. Do you know the cause of this increasing tread? As I post above, the build in IFFT doesn't have it. unfortunately, I have to use analytical equation to rebuild the signal from FFT components.       

Interesting... Have you tried to add the sin term instead of only cos?

Simplifying:

Amplitude * (cos(2*pi*freq*x+phase) + sin(2*pi*freq*x+phase))

samecues, thanks for your reply. I tried cos and exp(both has the same increasing tread), but not both sin and cos. I read there is discussion on internet "Using both sin and cosine terms will, unless there is something odd about the FFT, cause scaling issues and a phase shift"