Filtering the two frequencies with highest amplitudes of a signal in the frequency domain-CodePudding

I have tried to filter the two frequencies which have the highest amplitudes. Unfortunately the result is wrong (zero). Where is my error?

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.fftpack import ifft

data = np.loadtxt("profil.txt")
t = data[:,0]
x = data[:,1]

amplitudes = fft(x)
c = np.abs(amplitudes)
c_half = c[0:len(amplitudes)//2]
c_sorted = c_half.argsort()
i_2max = c_sorted[-2:]
p = np.zeros(len(amplitudes))
p[i_2max] = 0
p[-i_2max] = 0

x_filtered = ifft(amplitudes*p)

plt.plot(t,x_filtered, color='red')
plt.plot(t,x)
plt.show()

The resulting plot: enter image description here

CodePudding user response：

This part is very wasteful:

a = fft(x).real
b = fft(x).imag

You’re computing the FFT twice for no good reason. You compute it a 3rd time later, and you already computed it once before. You should compute it only once, not 4 times. The FFT is the most expensive part of your code.

Then:

ab_filter_2 = fft(x)
ab_filter_2.real = a*p
ab_filter_2.imag = b*p
x_filter2 = ifft(ab_filter_2)*2

Replace all of that with:

out = ifft(fft(x) * p)

Here you do the same thing twice:

p[c[0:int(len(c)/2)].argsort()[int(len(c)/2-1)]] = 0
p[c[0:int(len(c)/2)].argsort()[int(len(c)/2-2)]] = 0

But you set only the left half of the filter. It is important to make a symmetric filter. There are two locations where abs(f) has the same value (up to rounding errors!), there is a positive and a negative frequency that go together. Those two locations should have the same filter value (actually complex conjugate, but you have a real-valued filter so the difference doesn’t matter in this case).

I’m unsure what that indexing does anyway. I would split out the statement into shorter parts on separate lines for readability.

I would do it this way:

import numpy as np

x = ...
x -= np.mean(x)
fft_x = np.fft.fft(x)
c = np.abs(fft_x)  # no point in normalizing, doesn't change the order when sorting

f = c[0:len(c)//2].argsort()
f = f[-2:] # the last two elements are the indices to the largest two frequency components
p = np.zeros(len(c))
p[f] = 1   # preserve the largest two components
p[-f] = 1  # set the same components counting from the end

out = np.fft.ifft(fft_x * p).real
# note that np.fft.ifft(fft_x * p).imag is approximately zero if the filter is created correctly

Is it correct that the output of the FFT-function contains the fundamental frequency A0/ C0 […]?

In principle yes, but you subtracted the mean from the signal, effectively setting the fundamental frequency (DC component) to 0.