Computing FFT of a spectrum using python

Question

The spectrum shows ripples that we can visually quantify as ~50 MHz ripples. I am looking for a method to calculate the frequency of these ripples other than by visual inspection of thousands of spectra. Since the function is in frequency domain, taking FFT would get it back into time domain (with time reversal if I am correct). How can we get frequency of these ripples?

Answer 1

The problem arises from the fact that you are making a confusion between the term 'frequency' which you are measuring and the frequency of your data.

What you want is the ripple frequency, which actually is the period of your data.

With that out of the way, let's have a look at how to fix your fft.

As pointed out by Dmitrii's answer, you must determine the sampling frequency of your data and also get rid of the low frequency components in your FFT result.

To determine the sampling frequency, you can determine the sampling period by subtracting each sample by its predecessor and computing the average. The average sampling frequency will just be the inverse of that.

fs = 1 / np.mean(freq[1:] - freq[:-1])

For the high pass filter, you may use a butterworth filter, this is a good implementation.

# Defining a high pass filter
def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y

Next, when plotting the fft, you need to take the absolute value of it, that is what you are after. Also, since it gives you both the positive and negative parts, you can just use the positive one. As far as the x-axis is concerned, it will be from 0 to half of your sampling frequency. This is further explored on this answer

fft_amp = np.abs(np.fft.fft(amp, amp.size))
fft_amp = fft_amp[0:fft_amp.size // 2]
fft_freq = np.linspace(0, fs / 2, fft_amp.size)

Now, to determine the ripple frequency, simply obtain the peak of the FFT. The value you are looking for (around 50MHz) will be the period of the ripple peak (in GHz), since your original data was in GHz. For this example, it is actually around 57MHz.

peak = fft_freq[np.argmax(fft_amp)]

ripple_period = 1 / peak * 1000

print(f'The ripple period is {ripple_period} MHz')

And here is the complete code, which also plots the data.

import numpy as np
import pylab as plt
from scipy import signal as signal


# Defining a high pass filter
def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y


with open('ripple.csv', 'r') as fil:
    data = np.genfromtxt(fil, delimiter=',', skip_header=True)

amp = data[:, 0]
freq = data[:, 1]


# Determine the sampling frequency of the data (it is around 500 Hz)
fs = 1 / np.mean(freq[1:] - freq[:-1])

# Apply a median filter to remove the noise
amp = signal.medfilt(amp)

# Apply a highpass filter to remove the low frequency components 5 Hz was chosen
# as the cutoff fequency by visual inspection. Depending on the problem, you
# might want to choose a different value

cutoff_freq = 5
amp = butter_highpass_filter(amp, cutoff_freq, fs)

_, ax = plt.subplots(ncols=2, nrows=1)
ax[0].plot(freq, amp)
ax[0].set_xlabel('Frequency GHz')
ax[0].set_ylabel('Intensity dB')
ax[0].set_title('Filtered signal')

# The FFT part is as follows

fft_amp = np.abs(np.fft.fft(amp, amp.size))
fft_amp = fft_amp[0:fft_amp.size // 2]
fft_freq = np.linspace(0, fs / 2, fft_amp.size)

ax[1].plot(fft_freq, 2 / fft_amp.size * fft_amp, 'r-')  # the red plot
ax[1].set_xlabel('FFT frequency')
ax[1].set_ylabel('Intensity dB')

plt.show()

peak = fft_freq[np.argmax(fft_amp)]

ripple_period = 1 / peak * 1000

print(f'The ripple period is {ripple_period} MHz')

And here is the plot:

Answer 2

To get a proper spectrum for the blue plot you need to do two things:

1. Properly calculate frequencies for the spectrum plot (the red one)
2. Remove bias in the data so the spectrum is less contaminated with low frequencies. That's because you're interested in the ripple, not in the slow fluctuations.

Note, that when you compute fft, you get complex values that contain information about both amplitude and phase of oscillations for each frequency. In your case, the red plot should be an amplitude spectrum (compared to the phase spectrum). To get that, we take absolute values of fft coefficients.

Also, the spectrum you get with fft is two-sided and symmetric (since the signal is real). You really need only one side to get the idea where your ripple peak frequency is. I've implemented this in code.

After playing with your data, here's what I've got:

import pandas as pd
import numpy as np
import pylab as plt
import plotly.graph_objects as go
from scipy import signal as sig

df = pd.read_csv("ripple.csv")
f = df.Frequency.to_numpy()
data = df.Data
data = sig.medfilt(data)  # median filter to remove the spikes

fig = go.Figure()
fig.add_trace(go.Scatter(x=f, y=(data - data.mean())))
fig.update_layout(
    xaxis_title="Frequency in GHz", yaxis_title="dB"
)  # the blue plot with ripples
fig.show()

# Remove bias to get rid of low frequency peak
data_fft = np.fft.fft(data - data.mean())

L = len(data)  # number of samples

# Compute two-sided spectrum
tssp = abs(data_fft / L)

# Compute one-sided spectrum
ossp = tssp[0 : int(L / 2)]
ossp[1:-1] = 2 * ossp[1:-1]

delta_freq = f[1] - f[0]  # without this freqs computation is incorrect
freqs = np.fft.fftfreq(f.shape[-1], delta_freq)

# Use first half of freqs since spectrum is one-sided
plt.plot(freqs[: int(L / 2)], ossp, "r-")  # the red plot
plt.xlim([0, 50])
plt.xticks(np.arange(0, 50, 1))
plt.grid()
plt.xlabel("Oscillations per frequency")
plt.show()

So you can see there are two peaks: low-freq. oscillations between 1 and 2 Hz and your ripple at around 17 oscillations per GHz.