2

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

Hiddenguy
  • 547
  • 13
  • 33
  • 1
    You can't "erase" noise; it's simply there, especially when you use the raw FFT for spectral density estimation. Smoothing is usually the way to go, so what makes you think that Welch does not work for you? The lower plot looks exactly like what I would expect Welch to produce. – MB-F Dec 04 '17 at 11:59
  • Cause welch looks like this (see question). Do you have any ideas then what to do with my data? – Hiddenguy Dec 04 '17 at 13:15
  • 1
    `welch` has an argument `nperseg`. Higher values give you better frequency resolution and lower values give more noise reduction. Start with `signal.welch(dataset, fs=300, window='hamming', nperseg=256)` then try `nperseg=512` if the noise is still acceptable. I suggest you play with the values to get a feeling for them (also have a look at the `noverlap` argument). – MB-F Dec 04 '17 at 13:35
  • @kazemakase As I said previously it does not work properly. First look at my graphs (No. 1 and 2) they show a big peak at 60 Hz, which means it is an electricity line in the US). Using Welch I get it at 20 Hz. And what is more important, welch makes a line, and I need something like in graph 2, which is obviously not a line, but a tiny wavey line, which for me is desired. – Hiddenguy Dec 04 '17 at 13:41
  • 1
    The frequency axes do not match. Either your FFT plots are plainly wrong or your sampling rate is not 300 Hz. (I hope it's the latter, otherwise this might indicate a problem with your data.) Have you tried playing with `nperseg` as I suggested? This should allow you to get any amount of tiny-wavyness (noise) you desire. – MB-F Dec 04 '17 at 14:27
  • Yeah I used what you have proposed but it didn't work either. I have even changed fs to 900, but also nothing has changed. I've got 266336 samples in each lines, time of measurement is 300 seconds. So I guess we are trying the wrong method, or maybe my explanation is misleading. – Hiddenguy Dec 04 '17 at 14:46
  • 1
    I see.. well, `fs` is the sampling frequency. If you have 266336 samples over 300 seconds it should be 266336 / 300. If you have that many samples You can try very high values of `noverlap`. Maybe 4096 or more? – MB-F Dec 04 '17 at 14:54
  • @kazemakase Look at updated question – Hiddenguy Dec 04 '17 at 15:08
  • 1
    Code and images are easier than comments :) See my answer that shows the effect of the `nperseg` parameter. (I accidentally wrote noverlap in the comment above.) It should be possible to achieve the effect you are looking for (8192 with my example data), no? – MB-F Dec 04 '17 at 15:14

1 Answers1

4

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()
MB-F
  • 22,770
  • 4
  • 61
  • 116