Search code examples
signal-processingpython-3.7window-functionsnoise

Power spectrum, Welch method and window functions


I need to find power spectrum of a signal applying Welch method. Then using Chebyshev, Kaiser and Gauss windows, find the best window in terms of frequency separation.

I've wrote a function that produces noised sinusoidal signal and used the Welch method to this signal, also generated all of these three window functions.

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

def sig_noise(f,snr,n):
    fs=1000
    dt=1/fs
    tk=n*dt
    A=1
    B=10**(-snr/20)
    t=np.arange(0,tk,dt)
    y=(np.random.rand(n)-0.5)
    z=A*np.sin(2*np.pi*f*t)+B*y
    return z,t

[s,t]=sig_noise(140,-10,64)

f, Pxx_den = sg.welch(s,nperseg=64,noverlap=None) # Welch method, parameters are taken from book

plt.plot(f,Pxx_den)
plt.xlabel("Frequency")
plt.grid()

plt.figure(figsize=(10,10))

chebyshev=sg.chebwin=(64,40) # Chebyshev window

kaiser=sg.kaiser(64,0) # Kaiser window

gauss=sg.gaussian(64,np.std((f,Pxx_den),ddof=0)) # Gauss window

# plots for check

plt.subplot(221)
plt.grid()
plt.title("Gauss")
plt.plot(gauss)
plt.subplot(222)
plt.grid()
plt.title("Chebyshev")
plt.plot(chebyshev)
plt.subplot(223)
plt.grid()
plt.title("Kaiser")
plt.plot(kaiser)

f, Pxx_den = sg.welch(s,nperseg=64,noverlap=None,window=chebyshev)

the last line of code doesn't work as I want to use for example chebyshev window. Unknown window type is what I get as the answer. Could somebody tell me where is the mistake?


Solution

  • The issue is with a extra = character on the following line:

    chebyshev=sg.chebwin=(64,40) # Chebyshev window
    #                  ^^^
    #                  extra "=" character
    

    This causes chebyshev to be equal to the tuple (64,40) rather than the intended array (following a similar approach as for your other windows). Later when given a tuple as windows argument, scipy.signal.welch tries to construct the window by calling scipy.signal.get_window for a window type specified by the first element of the tuple. Obviously 64 (the first element of the tuple in your case) is not a valid window type.

    So, the fix is to simply correct the typo:

    chebyshev=sg.chebwin(64,40) # Chebyshev window
    

    or alternatively call welch with an appropriately constructed tuple:

    f, Pxx_den = sg.welch(s,nperseg=64,noverlap=None,window=('chebwin',40)