The Situation

I am currently writing a program that will later on be used to analyze a signal that is somewhat of a asymmetric Gaussian. I am interested in how many frequencies I need to reproduce the signal somewhat exact and especially the amplitudes of those frequencies. Before I input the real data I'm testing the program with a default (asymmetric) Gaussian, as can be seen in the code below.

My Problem

To ensure I that get the amplitudes right, I am rebuilding the original signal using the whole frequency spectrum, but there are some difficulties. I get to reproduce the signal somewhat well multiplying amp with 0.16, which I got by looking at the fraction rebuild/original. Of course, this is really unsatisfying and can't be the correct solution.

To be precise the difference is not dependant on the time length and seems to be a Gaussian too, following the form of the original, increasing in asymmetry according to the Skewnorm function itself. The amplitude of the difference function is correlated linear to 'height'.

My Question

I am writing this post because I am out of ideas for getting the amplitude right. Maybe anyone has had the same / a similar problem and can share their solution for this / give a hint.

Further information

Before focusing on a (asymmetric) Gaussian I analyzed periodic signals and rectangular pulses, which sadly were very unstable to variations in the time length of the input signal. In this context, I experimented with window functions, which seemed to speed up the process and increase the stability, the reason being that I had to integrate the peaks. Working with the Gaussian I got told to take each peak, received via the bare fft and ditch the integration approach, therefore my incertitude considering the amplitude described above. Maybe anyone got an opinion on the approach chosen by me and if necessary can deliver an improvement.

Code

from numpy.fft import fft, fftfreq
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skewnorm
np.random.seed(1234)


def data():
    height = 1
    data = height * skewnorm.pdf(t, a=0, loc=t[int(N/2)])
    # noise_power = 1E-6
    # noise = np.random.normal(scale=np.sqrt(noise_power), size=t.shape)
    # data += noise

    return data


def fft_own(data):
    freq = fftfreq(N, dt)
    data_fft = fft(data) * np.pi
    amp = 2/N * np.abs(data_fft)  # * factor (depending on t1)
    # amp = 2/T * np.abs(data_fft)**2
    phase = np.angle(data_fft)
    peaks, = np.where(amp >= 0)  # use whole spectrum for rebuild

    return freq, amp, phase, peaks


def rebuild(fft_own):
    freq, amp, phase, peaks = fft_own

    df = freq[1] - freq[0]
    data_rebuild = 0
    for i in peaks:
        amplitude = amp[i] * df
        # amplitude = amp[i] * 0.1
        # amplitude = np.sqrt(amp[i] * df)
        data_rebuild += amplitude * np.exp(0+1j * (2*np.pi * freq[i] * t
                                                   + phase[i]))

    f, ax = plt.subplots(1, 1)
    # mask = (t >= 0) & (t <= t1-1)
    ax.plot(t, data_init, label="initial signal")
    ax.plot(t, np.real(data_rebuild), label="rebuild")
    # ax.plot(t[mask], (data_init - np.real(data_rebuild))[mask], label="diff")
    ax.set_xlim(0, t1-1)
    ax.legend()


t0 = 0
t1 = 10  # diff(t0, t1) ∝ df
# T = t1- t0
N = 4096
t = np.linspace(t0, t1, int(N))
dt = (t1 - t0) / N

data_init = data()
fft_init = fft_own(data_init)
rebuild_init = rebuild(fft_init)

You should get a perfect reconstruction if you divide amp by N, and remove all your other factors.

Currently you do:

data_fft = fft(data) * np.pi  # Multiply by pi
amp = 2/N * np.abs(data_fft)  # Multiply by 2/N
amplitude = amp[i] * df       # Multiply by df = 1/(dt*N) = 1/10

This means that you currently multiply by a total of pi * 2 / 10, or 0.628, that you shouldn't (only the 1/N factor in there is correct).

Correct code:

def fft_own(data):
    freq = fftfreq(N, dt)
    data_fft = fft(data)
    amp = np.abs(data_fft) / N
    phase = np.angle(data_fft)
    peaks, = np.where(amp >= 0)  # use whole spectrum for rebuild
    return freq, amp, phase, peaks

def rebuild(fft_own):
    freq, amp, phase, peaks = fft_own
    data_rebuild = 0
    for i in peaks:
        data_rebuild += amp[i] * np.exp(0+1j * (2*np.pi * freq[i] * t
                                                + phase[i]))

Your program can be significantly simplified by using ifft. Simply set to 0 those frequencies in data_fft that you don't want to include in the reconstruction, and apply ifft to it:

data_fft = fft(data)
data_fft[np.abs(data_fft) < threshold] = 0
rebuild = ifft(data_fft).real

Note that the Fourier transform of a Gaussian is a Gaussian, so you won't be picking out individual peaks, you are picking a compact range of frequencies that will always include 0. This is an ideal low-pass filter.