CODEWITHSUNDEEP
pythonnumpyscipyfft
rmiller415
rmiller415

Reputation: 43

Manually filtering a FFT by frequency returns a shifted plot

I am taking a numerical and computational class, so we are asked to do things that can be done more efficiently by other packages, but we have to do it how we are asked. In this problem I was given some code that makes a noisy gaussian and asked to filter above or below certain frequencies. So, I got the fft frequencies and used list comprehension to set all of the frequencies greater than the absolute value of 2. Then, I use a for loop to set the coefficients to zero when the frequencies are zero. Then perform an inverse fft to get a filtered gaussian back. I get the correct graph out, but it's shifted down by 0.03 and it seems like the for loop is doing it. I'm not sure why though. Here is the code and the figure I get out.

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


##############Start Original Code###############
def gaussNoisy(x, noiseAmp):
    noise=noiseAmp*(np.random.randn(len(x)))
    return np.exp(-0.5*x*x)*(1.0+noise)
N=256
x=np.linspace(-4.0, 4.0, N)
y=gaussNoisy(x, 0.1)
##################End Original Code############

dx = x[1] - x[0]#define step size
Y_x = fft.fft(y) #perform fft
Y_x_shift = fft.fftshift(Y_x)
FFT_freq = fft.fftfreq(y.size,d=dx) #frequency spectrum
FFT_freq_shift = fft.fftshift(FFT_freq)

fig = plt.figure(1, figsize=(16,12))
ax1 = fig.add_subplot(211)
ax1.plot(x,y)
ax1.set_xlabel('Position')
ax1.set_ylabel('Amplitude')







#FFT_filtered32 = np.copy(FFT_freq) #Initialize FFT filter arrays
#FFT_filtered8 = np.copy(FFT_freq)
FFT_filtered2 = np.copy(FFT_freq)

Y_x2 = np.copy(Y_x)
Y_x32 = np.copy(Y_x)



FFT_filtered32 = [item if -32 <= item <= 32 else 0 for item in FFT_filtered32] #Filters arrays between +/- 32, 8, 2
FFT_filtered8 = [item if -8 <= item <= 8 else 0 for item in FFT_filtered8]
FFT_filtered2 = [item if -2 <= item <= 2 else 0 for item in FFT_filtered2]




l = 0
for l in range(Y_x.size):
    if  FFT_filtered2[l] == 0:
        Y_x2[l] = 0
    l += 1
    

fig = plt.figure(1, figsize=(16,12))
ax1 = fig.add_subplot(211)
ax1.plot(x,y)
ax1.set_xlabel('Position')
ax1.set_ylabel('Amplitude')
ax1.plot(x, fft.ifft(Y_x2))

Original noisy (blue) and filtered (orange).

Upvotes: 2

Views: 288

Answers (1)

Pascal Getreuer
Pascal Getreuer

Reputation: 3256

The cause of this problem is that the DC component is getting zeroed out in the l loop. This causes the apparent 0.3 shift of the signal, since the original is not zero mean. To fix this, skip l=0 for instance by starting the loop at 1: l in range(1, Y_x.size).

Upvotes: 1

Related Questions