Spectral estimation

Spectral estimation#

Given a time series

(6)#\[\begin{equation} x(n) \end{equation}\]

with fft spectrum

(7)#\[\begin{equation} y(m) = \sum_{n=0}^{N-1}x(n) \exp\{-2\pi i\frac{nm}{N}\} \end{equation}\]

or

(8)#\[\begin{equation} y(m) = FFT(x(n)) \end{equation}\]

Power spectrum is defined as

(9)#\[\begin{equation} P(m) = |y(m)|^2 \end{equation}\]

Power spectral density (PSD) of a discrete-time noise signal is given by the FFT of its autocorrelation function, R(k)

(10)#\[\begin{equation} PSD(m) = FFT(R(k)) \end{equation}\]

where the autocorrelation function R(k) is defined as

(11)#\[\begin{equation} R(k)=\frac{1}{N}\sum_{n=0}^{N-1}x(n)x(n-k) \end{equation}\]

that is

(12)#\[\begin{equation} PSD(m) = \sum_{k=0}^{N-1}\frac{1}{N}\sum_{n=0}^{N-1}x(n)x(n-k) \exp\{-2\pi i\frac{km}{N}\} \end{equation}\]

(13)#\[\begin{equation} PSD(m) = \frac{1}{N}\sum_{n=0}^{N-1}x(n)\exp\{-2\pi i\frac{nm}{N} \} \sum_{k=0}^{N-1}x(n-k) \exp\{2\pi i\frac{(n-k)m}{N}\} \end{equation}\]

(14)#\[\begin{equation} PSD(m) = \frac{1}{N}|FFT(x(n))|^2 \end{equation}\]

or

(15)#\[\begin{equation} PSD(m) = \frac{1}{N} P(m) \end{equation}\]

While the PSD is related to the Power Spectrum by dividing by length of FFT (i.e. length of signal window), the proper estimation of the Autocorrelation funtion \(R(k)\) requires a signal length that is typically longer than the FFT, which cannot be done for short signals (transients), even if mathematically possible, but estimate is not reliable.

1/Hz density estimate#

Above PSD is on a per sample bases. Concerning the 1/Hz density estimate one has to scale the estimate by the number of samples per second.

(16)#\[\begin{equation} PSD(m) = \frac{1}{fs N} P(m) \end{equation}\]

where \(fs\) is the sampling frequency

Windowed PSD estimates#

When windowing the timeseries before PSD estimates

(17)#\[\begin{equation} PSD(m) = \frac{1}{fs \sum{w(n)^2}} |FFT(w(n)x(n))|^2 \end{equation}\]

For a rectangular window we have \(\sum{w(n)^2}=N\)

Windowed power spectrum estimates#

Alternatively to the PSD we have the power spectrum

(18)#\[\begin{equation} PS(m) = \frac{1}{(\sum{w(n)})^2} |FFT(w(n)x(n))|^2 \end{equation}\]

that allows the RMS spectral estimates of tonal signals, but provides wrong results for random signals

Welch PSD#

In the following, some Python should demonstrate the different

Startup and helper functions.#

There are two functions defined

mSpectrum: direct FFT based power spectral estimate
wSpectrum: Welch based estimation (for comparison as well defind method)
estSpectrum: spectrum estimation and plotting function

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

def mSpectrum(xx,win,nfft,fs):
    # estimate power spectral density (PSD) usinf rfft directly
    nwin=len(win)
    yy=np.fft.rfft(xx[:nwin]*win,nfft)
    pp=np.abs(yy)**2
    pp[1:] *=2                  # compensate for negatice frequencies
    ff=np.arange(int(nfft/2)+1)*fs/nfft
    sc1= fs * (win*win).sum()  # for density
    sc2= win.sum()**2          # for spectrum
    return ff,pp/sc1,pp/sc2

def wSpectrum(xx,win,nfft, fs):
    # estimate power spectral density (PSD)
    f, Pxx_dens = signal.welch(xx, fs=fs, nperseg=len(win), nfft=nfft, noverlap=0, scaling="density", window=win)
    # estimate power spectrum
    f, Pxx_spec = signal.welch(xx, fs=fs, nperseg=len(win), nfft=nfft, noverlap=0, scaling="spectrum", window=win)
    return f,Pxx_dens,Pxx_spec


def estSpectra(time,xx,win,nfft,fs,amp,xl,yl):
    f,Pxx_dens,Pxx_spec=mSpectrum(xx,win,nfft,fs)
    #
    print("'PSD' peak (RMS):      ",np.sqrt(Pxx_dens.max()), ';  Noise:', np.mean(Pxx_dens[3*len(f)//4:]))
    print("'Spectrum' peak (RMS): ",np.sqrt(Pxx_spec.max()), ';  Noise:', np.mean(Pxx_spec[3*len(f)//4:]))
    #
    fig = plt.figure("figure.figsize",[15,10])

    plt.subplot(211)
    plt.plot(time,xx)
    plt.plot(time,amp+0*time,'r')
    plt.xlim(xl)
    plt.grid(True)

    plt.subplot(223)
    plt.semilogy(f/1000, Pxx_dens)
    plt.ylim(yl)
    plt.xlabel('Frequency [kHz]')
    plt.ylabel('[V**2/Hz]')
    plt.title('PSD')
    plt.grid(True)

    plt.subplot(224)
    plt.semilogy(f/1000, Pxx_spec)
    plt.ylim(yl)
    plt.xlabel('Frequency [kHz]')
    plt.ylabel('[V**2]')
    plt.title('Spectrum')
    plt.grid(True)

Tonal + Noise#

The sampling frequency, signal and window length are chosen such that the Welch PSD corresponds to a single FFT transform. For longer signals the Welch PSD estimates would only average the spectra.

From the graphs one notes that the noise level (0.001) is estimated by the PSD (density) method, while the signal RMS (10 units) estimate corresponds to the peak value of the Spectrum (not density) estimate, but the noise density is then wrongly estimated.

The examples are derived from the Python Welch documentation and are consequently only valid in the python environment.

# simulate tone
def genSignal(freq,fs,N=None,arms=1,noise=0):
    # f: signal frequency
    # fs: sampling frequency
    # N: number of samples
    # arms: rms amplitude of signal
    # noise: noise variance
    if N==None: N=fs
    #
    time=np.arange(N)/fs
    #
    amp=arms*np.sqrt(2)   
    noise_power=noise *fs/2 # convert to noise bandwidth (-1/T to 1/T)

    xx = amp*np.sin(2*np.pi*freq*time)
    #rms= np.sqrt(np.mean(xx**2))
    xx += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
    return time,xx

fs = 1024   #Hz
N = fs*10

#    Generate a 1 second test signal, a 10 Vrms sine wave at fs/8, corrupted by
#    0.001 V**2/Hz of white noise sampled at  2.048 kHz.

freq = fs/8

arms=10
noise=0.001
time,xx=genSignal(freq,fs,N,arms,noise)

print("Signal RMS :",arms, ';  noise :', noise)

# process and display
nwin=N
nfft=nwin*2
# obtain window
win=signal.get_window("hann",nwin)
#
estSpectra(time,xx,win,nfft,fs,arms,[0.15,0.35],[1e-6,1e+3])

Signal RMS : 10 ;  noise : 0.001
'PSD' peak (RMS):       25.78624185281698 ;  Noise: 0.0009822312853819209
'Spectrum' peak (RMS):  9.986968525723686 ;  Noise: 0.00014733469280728815

_images/9cb5cb552785817692635c6da65801af19b48af98413032c72394ec4bea1aac8.png

Transient + noise#

Assume short transient (echolocation click) in the middle of the time window

The transient is sinusoidal with a bell shaped (gaussian) amplitude function.

We can see from the pictures that the noise density is correctly estimated by the PSD (density) method, but both the density and the spectrum method cannot provide a good description of the transient. The signal is to short for RMS estimation. Signal characteristics can only be estimated after signal detection, when the background noise is elimiated.

fs = 1024   #Hz

#    Generate a 1 second test signal, a 10 VPeak Gaussian shaped sine wave at fs/4, corrupted by
#    0.001 V**2/Hz of white noise sampled at  1.024 kHz.
T=2
N = fs*T
tt = np.arange(N) / fs
to = N//(2*T)
sigma=8/fs
freq = fs/4

arms=40*np.exp(-1/2*((tt-tt[to])/sigma)**2)
noise=0.001
time,xx=genSignal(freq,fs,N,arms,noise)

print("Signal peak: ",arms.max())#, energy)


nwin=N
nfft=nwin*2


#time = np.arange(N) / fs

#amp=10*np.exp(-1/2*((time-time[N//2])/sigma)**2)

#freq = fs/4
#noise_power = 0.001 * fs / 2 
#xx = amp*np.sin(2*np.pi*freq*time) 
#energy= np.sum(xx**2)
#xx += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

# obtain window
win=signal.get_window("hann",nwin)

estSpectra(time,xx,win,nfft,fs,arms,[0.4,0.6],[1e-6,1e+3])

Signal peak:  40.0
'PSD' peak (RMS):       0.4865374833512813 ;  Noise: 0.0011052008843724917
'Spectrum' peak (RMS):  0.42135382047555797 ;  Noise: 0.0008289006632793687

_images/b898beca6162df340422d6cd5f0496f71c698ed7399d4ca12d00c1809b7a2bf7.png

This examples shows that when window or FFT length is much longer than signal, the contamination of noise is significant and does not allow a reliable spectral description of the the signal.

Estimmation of transient peak amplitude#

After shown that both metrics V^2/Hz and RMS make no sense for tansients and therefore echolcation clicks, the question arizes how to scale the spectrum in order to obtain spectral peak values that correspond to the time series.

First we note that the signal is smaller than the analysis window so we need to remove contributions from outside the signal from the spectrum. This can done by either a rectangular that is result of the signal detection. In the following a gaussian window is used.

Assume signal with gaussian shape

(19)#\[\begin{equation} x(t)=\exp{(-at^2)} \end{equation}\]

Fourier transform

(20)#\[\begin{equation} X(\omega)=\int_{-\infty}^{\infty}{x(t) e^{-i\omega t} dt} \end{equation}\]

(21)#\[\begin{equation} X(\omega)=\int_{-\infty}^{\infty}{ \exp{-at^2-i\omega t} dt}=\int_{-\infty}^{\infty}{ \exp{-(at^2+i\omega t)} dt} \end{equation}\]

With

(22)#\[\begin{equation} (at^2+i\omega t)=at^2+2\sqrt{a} t \frac{i \omega}{2\sqrt{a}} +(\frac{i \omega}{2\sqrt{a}})^2-(\frac{i\omega}{2\sqrt{a}})^2 \end{equation}\]

or

(23)#\[\begin{equation} (at^2+i\omega t)=\Big(\sqrt{a}t+\frac{i \omega}{2\sqrt{a}}\Big)^2+(\frac{\omega}{2\sqrt{a}})^2 \end{equation}\]

(24)#\[\begin{equation} X(\omega)=\exp{(-\frac{\omega^2}{4a})}\int_{-\infty}^{\infty}{ \exp{-\Big(\sqrt{a}t+(\frac{i\omega}{2\sqrt{a}} )\Big)^2} dt} \end{equation}\]

with

(25)#\[\begin{equation} u(t)= \sqrt{a}t+(\frac{i\omega}{2\sqrt{a}} ) \end{equation}\]

and

(26)#\[\begin{equation} dt=\frac{du}{\sqrt{a}} \end{equation}\]

(27)#\[\begin{equation} X(\omega)=\exp{(-\frac{\omega^2}{4a})}\int_{-\infty}^{\infty}{ \exp{\big[-u(t)^2]} \frac{du}{\sqrt{a}}} = \sqrt{\frac{\pi}{a}}\exp{(-\frac{\omega^2}{4a})} \end{equation}\]

Gaussian modulated signal#

(28)#\[\begin{equation} x(t)=\exp{(-at^2)}\exp{(i\omega_0 t)} \end{equation}\]

using frequency shifting property of fourirt transforms:

if

(29)#\[\begin{equation} \text{FFT}\{x(t)\}=X(\omega) \end{equation}\]

then

(30)#\[\begin{equation} \text{FFT}\{x(t)\exp(i\omega_0 t)\}=X(\omega-\omega_0) \end{equation}\]

Spectral peak#

spectral peak of a gaussian amplitude is with

(31)#\[\begin{equation} a=\frac{1}{2\sigma^2} \end{equation}\]

(32)#\[\begin{equation} \sqrt{\frac{\pi}{a}}=\sqrt{2\pi}\sigma \end{equation}\]

fs=2000
T=5
N = fs*T
tt = np.arange(N) / fs
to = N//(2*T)
sigma=8/fs
freq = fs/4

nwin=N
aa=100
ashape=np.exp(-1/2*((tt-tt[to])/sigma)**2)
#ashape=signal.get_window("hann",nwin)

arms=aa*np.sqrt(2)*ashape
time,xx=genSignal(freq,fs,N,arms,0)
#
nfft=2*nwin
fr=np.arange(1+nfft//2)*fs/nfft
P =np.abs(np.fft.rfft(xx,nfft)/fs)
Pn=np.abs(np.fft.rfft(ashape,nfft)/fs).max()
P=P/Pn
print(P.max())

plt.plot(fr,P)
plt.grid(True)
plt.show()

100.00000000000027

_images/af9b4f5352dfda08d18e20fe96c13e1c83a4f04cabd609eaa0d78402f7c1a987.png

We see that scaling the power spectrum by the peak value of the spectrum of the pulse amplitude function normalizes the click spectrum so that spectral peak corrsponds to the time series peak anplitude.

Summary#

the above examples show that for different type of signals one uses different characteristic measurements resulting in different algorithms

noise spectral density#

Scale FFT power estimate with \(1/(f_s \sum w^2[n])\) to obtain power/Hz

tonals RMS estimates#

Scale FFT peak power estimate with \(1/(\sum w[n])^2\) to obtain signal rms estimate

transient peak estimates#

Scale FFT peak power estimate with peak of the spectrum of the transient amplitude function