Search code examples
pythonpython-2.7numpyfftnoise

Erasing noise from fft chart

Do you know how to delete so much noise from the FFT? Here is my code of FFT:

import numpy as np

fft1 = (Bx[51:-14])
fft2 = (By[1:-14])


# Loop for FFT data
for dataset in [fft1]:
    dataset = np.asarray(dataset)
    psd = np.abs(np.fft.fft(dataset))**2
    freq = np.fft.fftfreq(dataset.size, float(300)/dataset.size)
    plt.semilogy(freq[freq>0], psd[freq>0]/dataset.size**2, color='r')


for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    psd2 = np.abs(np.fft.fft(dataset2))**2
    freq2 = np.fft.fftfreq(dataset2.size, float(300)/dataset2.size)
    plt.semilogy(freq2[freq2>0], psd2[freq2>0]/dataset2.size**2, color='b')

What I get: enter image description here

What I need: enter image description here

Any ideas? Welch does not work, so as you can see, I don't want to smooth my chart, but erase so much noise to the level which is presented on the second picture.

This is what Welch do: enter image description here and a bit of code:

freqs, psd = scipy.signal.welch(dataset, fs=300, window='hamming')

Updated Welch: enter image description here

A bit of code:

# Loop for FFT data
for dataset in [fft1]:
    dataset = np.asarray(dataset)
    freqs, psd = welch(dataset, fs=266336/300, window='hamming', nperseg=512)
    plt.semilogy(freqs, psd/dataset.size**2, color='r')

for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    freqs2, psd2 = welch(dataset2, fs=266336/300, window='hamming', nperseg=512)
    plt.semilogy(freqs2, psd2/dataset2.size**2, color='b')

As you can see Welch is well configurated, it shows 60 Hz electricity line, and harmonic modes. It is almost good, but it smoothed completely my plot. See graph two which is desired. Btw. y scale is wrong at Welch plot, but it is just a case of power data to the two.

I have changed to nperseg=8192 and it worked. Look at the results. enter image description here

Solution

  • Here is an example that shows how to use nperseg to control the frequency resolution vs. noise reduction tradeoff:

    enter image description here

    Setting nperseg to the length of the signal is more or less equivalent to using the FFT without any averaging.

    Here is the code to generate this image:

    import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

plt.figure(figsize=[8, 12])

n = 2**21
fs = 887

# example data
x = np.random.randn(n)
x += np.sin(np.cumsum(0.42 + np.random.randn(n) * 0.01)) * 5
x = signal.lfilter([1, 0.5], 2, x)

plt.subplot(3, 2, 1)
plt.semilogy(np.abs(np.fft.fft(x)[:n//2])**2 / n**2, label='FFT')
plt.legend(loc='best')

for i, nperseg in enumerate([128, 512, 8192, 65536, n]):
    plt.subplot(3, 2, i+2)
    f, psd = signal.welch(x, fs=fs, window='hamming', nperseg=nperseg, noverlap=0)
    plt.semilogy(f, psd, label='nperseg={}'.format(nperseg))
    plt.legend(loc='best')

plt.show()