Hamming Window, python 2.7

Question

Hi I have a FFT which is quite noisy. How to apply to my code Hamming window to make it less noisy. Look at my code:

plt.subplot(212)
plt.title('Fast Fourier Transform')
plt.ylabel('Power [a.u.]')
plt.xlabel('Frequency Hz')
fft1 = (Bx[51:-14])
fft2 = (By[1:-14])

for dataset in [fft1]:
    dataset = np.asarray(dataset)
    psd = np.abs(np.fft.fft(dataset))**2.5
    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.5
    freq2 = np.fft.fftfreq(dataset2.size, float(300)/dataset2.size)
    plt.semilogy(freq2[freq2>0], psd2[freq2>0]/dataset2.size**2, color='b')

What plt.show() is

What I need is:

I have seen (https://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.signal.hamming.html) and this (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.hamming.html) but still don't have a clue how to apply it to my code. Any ideas? As I said you see in the second picture what I need. Maybe Blackman window will also be good to apply, but still not a clue how to add it.

Applying this:

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

Gave me that, which does not appear like my desired chart.

It looks like 300 isn't the actual sample rate. I was guessing based on your question, so use whatever is your actual sampling rate (could be around 900?). Also, keep in mind that Welch's method computes a smoother estimate of the power spectrum by sacrificing frequency resolution. — bnaecker, Dec 04 '17 at 04:08

bnaecker · Answer 1 · 2017-12-03T20:53:37.283

4

It seems like you're trying to estimate the power spectrum of your signals. If that's the case, you can use something like the scipy.signal.welch function, which estimates a smoothed spectrum by computing FFTs from overlapping segments of the data. You can pass the method a window keyword-argument directly, such as 'hamming' or 'blackman'.

EDIT:

Applying this to your data, you'd do something like this:

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

This will return the frequencies and power at those frequencies. I'm assuming here that 300 is your sample rate (from your computation of freq in your question).

edited Dec 03 '17 at 20:53

answered Dec 03 '17 at 15:48

bnaecker

6,152
1
20
33

@Hiddenguy What do you mean? You can look at the documentation for details on how to use it. Briefly, pass it a signal and the sampling rate, and it will return the frequencies and power at those frequencies. – bnaecker Dec 03 '17 at 17:15
I applied to my question two webpages of the documentation, but I have no idea how to implement them to my code. So can you help? – Hiddenguy Dec 03 '17 at 17:41
Using Welch is pointless, because the results always going to be wrong. I think my method is correct, all I need is to erase that much noise from my chart. – Hiddenguy Dec 04 '17 at 09:36
@Hiddenguy The results are not "always going to be wrong." Welch's method was literally *created* for the problem you're trying to solve. You may not be using the correct sampling rate, and you may need to play with the various keyword parameters to get the results you want. You can also try directly smoothing the power spectrum you compute in your question, but that is going to be sub-optimal compared to Welch's method. If you go that route, you're not computing the power correctly. It should be `np.square(np.abs(np.fft.fft(signal)))`. You're using 2.5 as the exponent, it should be 2.0. – bnaecker Dec 04 '17 at 16:47

score 0 · Accepted Answer · answered Dec 30 '17 at 19:07

It seems that welch method was correct so I thought about my question and here is the answer for it.

   # Loop for FFT data
    for dataset in [fft1]:
        dataset = np.asarray(dataset)
        freqs, psd = welch(dataset, fs=266336/300, window='hamming', nperseg=8192)
        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=8192)
        plt.semilogy(freqs2, psd2/dataset2.size**2, color='b')

Hamming Window, python 2.7

2 Answers2