BLM_1: Basic signal spectral analysis

BLM_1: Basic signal spectral analysis#

The following shows some basic steps that one undertakes for signal processing.

# loading basic modules
import numpy as np
import matplotlib.pyplot as plt

Simulation of sound wave

# creation of time vector
fs = 1024
nd = int(10*fs)
tt = np.arange(nd)/fs

# creation of cosine signal
amp=2
xx = amp*np.cos(2*np.pi*tt)

# plot signal
fig = plt.figure("figsize",[15,5])
plt.plot(tt,xx);
plt.show();

../_images/006b1cbef06e1f56bbdc614771e8960c3c88ef789e6b827f35ff6165d40ff5ec.png

#creation of sine signal
yy = amp*np.sin(2*np.pi*tt)

fig = plt.figure("figsize",[15,5])
plt.plot(tt,xx);
plt.plot(tt,yy);
plt.show();

../_images/be81e5256e0115ee9375821875bbc9e58a8ba4b56b266116a86300c60c9dccc9.png

#estimation of mangnitude
zz=np.sqrt(xx**2+yy**2)

fig = plt.figure("figsize",[15,5])
plt.plot(tt,xx);
plt.plot(tt,yy);
plt.plot(tt,zz);
plt.show();

../_images/9503c0a6b2b284b7db28434244858abe8ebae691337f6b639d2db04aab2310a5.png

print(np.mean(xx),np.mean(yy))
print(np.mean(xx**2),np.mean(yy**2))

-8.326672684688674e-17 -1.1102230246251566e-17
1.9999999999999996 2.0

Observations#

the mean of a sinosoidal wave is zero
the mean of the signals squared is equal to the amplitude.

Complex numbers#

The following mathematical excurse is helpful for further spectral analysis

Trigonometry#

Triangle in circle

import numpy as np
import matplotlib.pyplot as plt
arg=np.arange(360)*np.pi/180
plt.plot(np.cos(arg),np.sin(arg),'--')
plt.plot([0,np.cos(np.pi/6)],[0,0],'k')
plt.plot([np.cos(np.pi/6),np.cos(np.pi/6)],[0,np.sin(np.pi/6)],'k')
plt.plot([0,np.cos(np.pi/6)],[0,np.sin(np.pi/6)],'k')
plt.grid(True)
plt.xlim(-1.1,1.5)
plt.text(1.0,0.15,"$sin(phi)$")
plt.text(0.25,-0.15,"$cos(\phi)$")
plt.text(-0.25,0.4,"$cos^2(\phi)+sin^2(\phi)=1$")
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

<>:11: SyntaxWarning: invalid escape sequence '\p'
<>:12: SyntaxWarning: invalid escape sequence '\p'
<>:11: SyntaxWarning: invalid escape sequence '\p'
<>:12: SyntaxWarning: invalid escape sequence '\p'
C:\Users\zimme\AppData\Local\Temp\ipykernel_33888\1296300766.py:11: SyntaxWarning: invalid escape sequence '\p'
  plt.text(0.25,-0.15,"$cos(\phi)$")
C:\Users\zimme\AppData\Local\Temp\ipykernel_33888\1296300766.py:12: SyntaxWarning: invalid escape sequence '\p'
  plt.text(-0.25,0.4,"$cos^2(\phi)+sin^2(\phi)=1$")

../_images/285899e094a042c1cad0ba0fea179cd0bf4575bc035c06ac924c6f38e1b4b68a.png

Triangle

(195)#\[\begin{equation} x=\cos(\phi) \end{equation}\]

and

(196)#\[\begin{equation} y=\sin(\phi) \end{equation}\]

Pythagoras#

(197)#\[\begin{equation} x^2 + y^2 =r^2 \end{equation}\]

applied here

(198)#\[\begin{equation} \cos^2(\phi) + \sin^2(\phi) =1 \end{equation}\]

Problem#

Find a formula, so that sum of squares becomes constant (e.g. \(r^2\) or 1)

Hint#

We know \((a+b)*(a-b) = a^2 - b^2\), and more general \((a+c*b)*(a-c*b) = a^2 - c^2b^2\).

Solution#

Define a new symbol \(c^2 =i^2= -1\)

where \(i\) is defined such that

(199)#\[\begin{equation} i^2 = -1 \end{equation}\]

This way we get

(200)#\[\begin{equation} (\cos(\phi)+i\sin(\phi))(\cos(\phi)-i\sin(\phi))= \cos^2(\phi) + \sin^2(\phi) =1 \end{equation}\]

Definitions#

complex number (real + i*imaginary)#

(201)#\[\begin{equation} z=(\cos(\phi)+i\sin(\phi)) = \exp(i\phi)= e^{i\phi} \end{equation}\]

conjugate of a complex number#

(202)#\[\begin{equation} z^*=(\cos(\phi)-i\sin(\phi)) = \exp(-i\phi) = e^{-i\phi} \end{equation}\]

absolute squared complex value by signal multiplied with its conjugate#

(203)#\[\begin{equation} |z|^2 = z z^* = (\cos(\phi)+i\sin(\phi))(\cos(\phi)-i\sin(\phi)) \end{equation}\]

Alternative method#

(204)#\[\begin{equation} |z|^2 = e^{i\phi}e^{-i\phi} = e^{i\phi-i\phi} =e^0= 1 \end{equation}\]

Complex notation simplifies mathematics, but for physical phenomena only the real part is of relevance

Also from

(205)#\[\begin{equation} z z^* = 1 \end{equation}\]

follows

(206)#\[\begin{equation} z^* = \frac{1}{z} \end{equation}\]

A more general expression is

(207)#\[\begin{equation} z z^* = |z|^2 \end{equation}\]

so that

(208)#\[\begin{equation} z^* = \frac{|z|^2}{z} \end{equation}\]

or equivalently

(209)#\[\begin{equation} z^{-1} = \frac{z^*}{|z|^2} \end{equation}\]

# creation of complex signal (note 1j )
xx = amp*np.exp(-2*np.pi*1j*tt)

# plot signal
fig = plt.figure("figsize",[15,5])
plt.plot(tt,np.real(xx));   # cos(xx)
plt.plot(tt,np.imag(xx));   # sin(xx)
plt.plot(tt,np.abs(xx));    # sqrt(z * conjugate(z))
plt.show();

../_images/5ed33b143c710302567152c89ed42d7f4b9267587f2457515d12d9838faeaa70.png

General frequencies#

# creation of a 1 second time vector 
fs = 1024 # sampling per second
nd = int(1*fs)
tt = np.arange(nd)/fs
amp=3

# signal frequency cycles/second = Hz
fo = 5
phi=2*np.pi*fo*tt
xx=amp*np.exp(- 1j *phi)

# plot signal
fig = plt.figure("figsize",[15,5])
plt.plot(tt,np.real(xx));
plt.grid(True)
plt.show()

../_images/b706de9ca5b6ca262554aacebab8cbe871e544240fd72e1356a28669a09c78f0.png

x5=amp*np.exp(- 1j * 2*np.pi*5*tt)
w5=np.exp( 1j * 2*np.pi*5*tt)
y5=x5*w5
plt.plot(np.real(y5));

../_images/38977505c41d8e4869c1fbb30da055f965bb48267bb25a3f8a3ed3e7f7f1a049.png

amp=np.exp(-0.5*((tt-0.5)/0.1)**2)
plt.plot(amp)
x5=amp*np.exp(- 1j * 2*np.pi*5*tt)
w5=np.exp( 1j * 2*np.pi*5*tt)
y5=x5*w5
plt.plot(np.real(x5));
plt.plot(np.abs(y5));

../_images/5baf7f0c814404c0a8ece1277fe5d220a4166a8554738f9a0e1a0e806689ce17.png

amp=1
fo=10
plt.plot(amp)
x5=amp*np.exp(- 1j * 2*np.pi*fo*tt)
plt.plot(np.real(x5));
plt.show()
#
plt.figure('figsize',[15,5])

w4=np.exp( 1j * 2*np.pi*(fo-1)*tt)
y4=x5*w4
plt.subplot(131)
plt.plot(np.real(y4));
plt.plot(np.imag(y4));

w5=np.exp( 1j * 2*np.pi*fo*tt)
y5=x5*w5
plt.subplot(132)
plt.plot(np.real(y5));
plt.plot(np.imag(y5));

w6=np.exp( 1j * 2*np.pi*(fo+1)*tt)
y6=x5*w6
plt.subplot(133)
plt.plot(np.real(y6));
plt.plot(np.imag(y6));
plt.show()

../_images/e82461db329dbf9265291c4191114bab290d9689f30267d28b1195ed706dccf4.png

../_images/286c3f28d552a1a6e187f9a4d716546be764ded2dc2f1a40f8518866a8b787c0.png

Frequency analysis#

Discrete Fourier transform (DFT)#

xx=np.exp(- 1j * 2*np.pi*5*tt)

w4=np.exp(2*np.pi*1j *4*tt)
w5=np.exp(2*np.pi*1j *5*tt)
w6=np.exp(2*np.pi*1j *6*tt)

y4=np.sum(xx*w4)
y5=np.sum(xx*w5)
y6=np.sum(xx*w6)

print(y4,y5,y6)

(-2.842170943040401e-14+0j) (1024+2.958665228503174e-16j) (-8.526512829121202e-14+5.684341886080802e-14j)

xx=np.exp(- 1j * 2*np.pi*5*tt)

# matrix of spectral analysis vectors
W=np.exp(2*np.pi*1j *np.outer(np.arange(fs),tt))

yy1=W @ xx
#yy1=np.dot(W,xx)
#yy1=W.dot(xx)
plt.plot(np.abs(yy1),'.-');
plt.show()

yy2=W @ np.real(xx)
plt.plot(np.abs(yy2),'.-');
plt.show()

../_images/b7aa583ee110f9ec5f67951085a205b2fa941358eba6ab6bc196b35167df80dd.png

../_images/f08a2a4e4f6abf5b8b3376aab68eddec34d8288145bc80736e3b22fd52a2e65c.png

amp=np.exp(-0.5*((tt-0.5)/0.05)**2)
fo=50
xx=amp*np.exp(- 1j * 2*np.pi*fo*tt)

# matrix of spectral analysis vectors
ff=np.arange(fs)

W=np.exp(2*np.pi*1j *np.outer(ff,tt))

yy1=W @ xx
#yy1=np.dot(W,xx)
#yy1=W.dot(xx)
plt.plot(ff,np.abs(yy1),'.-');
plt.xlim(0,100)
plt.show()

yy2=W @ np.real(xx)
plt.plot(ff,np.abs(yy2),'.-');
plt.xlim(0,100)
plt.show()

../_images/53b2f6831150009789763e67a74e72f0c2299d145f4d04f3a658460b68e2eee2.png

../_images/bb872ddedf4ee590f5c4b1342bb89f333c705d422fd46366b471002a1c9a33d8.png

print(np.where(np.abs(yy1)>np.max(np.abs(yy1))/2))
print(np.where(np.abs(yy2)>np.max(np.abs(yy2))/2))

(array([47, 48, 49, 50, 51, 52, 53], dtype=int64),)
(array([ 47,  48,  49,  50,  51,  52,  53, 971, 972, 973, 974, 975, 976,
       977], dtype=int64),)

Observation#

complex signal (cos and sin) peak is at correct frequency (5) with peak amplitude that correspond to the length of the DFT
for real signals there are two peaks with half the amplitude at frequencies 5 and 1019 = (1024-5). The second peak is allso called negative frequency
the peak amplitude value of the spectrum of the complex signal is amp*length(fft)
the peak amplitude values of the spectrum of the real signal is (amp/2)*length(fft)

Fast Fourier Transform#

# generic complex to complex fft
yy=np.fft.fft(xx)

plt.plot(np.abs(yy),'.-');
plt.show()

# special real to complex fft
yy0=np.fft.rfft(np.real(xx))

plt.plot(np.abs(yy0),'.-');
plt.show()

../_images/e14313efae5308fc652e0f5ce9dc76f8c6d9d8ac37198ca1185979ac920b0d73.png

../_images/e1e11becc8cdc16557e01f7f6baca56da96e6801b6b73fef26754976970ad02d.png

Comments#

For a complex sinusoidal signal the Python FFT gives a negative frequency (fs-5 = 1019)
a special version of the FFT call rfft gives only the positive frequencies (x-axis only up to fs/2 and not fs)

FFT implementation#

To see how FFT is implemented (\(\exp(2\pi i nm/N)\) or \(\exp(-2\pi i nm/N)\)) we fourier transform a time series where only a single sample is set to a non-zero value (say sample 1 set to one)

In the case of a DFT we note that for a single sample

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

would become

(211)#\[\begin{equation} y(m)= \exp(2\pi i m/N)= \cos(2\pi m/N)+1i\sin(2\pi m/N) \end{equation}\]

taking the inverse DFT

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

which results to

(213)#\[\begin{equation} x(n) = M\delta_{n,1} \end{equation}\]

or \(x(1)=M\) and \(x(n)=0\) for \(n \ne 1\)

# 
v1 = np.exp(2*np.pi*1j*np.arange(1024)/1024)
plt.plot(np.real(v1))
plt.plot(np.imag(v1))
plt.show()

u1=np.sum(v1*np.exp(-2*np.pi*1j*np.arange(1024)/1024))
print("'DFT' peak: ", np.real(u1))

../_images/348a749ce99a12a5a9e5958b8a579e2f1be35ac6cc14ad18c98012ce126a8c61.png

'DFT' peak:  1024.0

uu=np.zeros(1024,complex)
uu[1]=1
vv=np.fft.fft(uu)

plt.plot(np.real(vv));
plt.plot(np.imag(vv));
plt.xlabel("Frequency");
plt.show()

ww=np.fft.ifft(vv)
print("FFT peak: ",np.real(uu[1]),np.real(ww[1]))

../_images/cb68d297c3418db8b710fab6fdbcce0eb3911b86b3444fac75918d4a6196afd9.png

FFT peak:  1.0 1.0

Observation#

The FFT spectrum of a signle sample is \(\cos(2\pi m/N) -i \sin(2\pi m/N)\) contrary to the DFT used above that was using \(\cos(2\pi m/N) +i \sin(2\pi m/N)\) . This selection \(+i\) or \(-1\) is purely a matter of cenvention. In practical terms it means that a decision has to be made if the ‘non-real’ imaginary signal is trailing (\(+i\)) or leading \(-i\) the real signal.

As we are dealing mostly real signals and use the special rfft that only results in positive frequencies, we have not to worry about this ambiguity.

While the spectra of the FFT and DFT are equivalent (equal when ignoring the different sign of the imaginary part), the inverse FFT produces the correct original time series, while the DFT

General signal frequency#

So far, the frequency was chosen to be an integer value. In the following we relax this and use a non-integer frequency, e.g.: \(f_o = 5.5\)

# signal frequency 
fo = 5.5
xx=amp*np.exp(-2*np.pi*1j *fo*tt)

# plot signal
fig = plt.figure("figsize",[15,5])
plt.plot(tt,np.real(xx));
plt.grid(True)
plt.show()

yy3=np.fft.rfft(np.real(xx))

fig = plt.figure("figsize",[15,5])
plt.plot(np.abs(yy3),'.-');
plt.show()

../_images/6f1b59ee60f836d547de3f6852365a7c7f374a5dd6fe621ebe93624a9271fd9b.png

../_images/f4197a840022b84addfcd6833d37528020f715995a166c91f775efee2a52c35a.png

Observations#

Using a fractional frequency (5.5 instead of 5) different effects occur

peak is reduced to less than 700
peak is wider
most frequencies are not zero

Conclusions#

When using fractional frequencies the energy (power) of the signal is leaked to other frequencies.

Spectral leakeage#

The widening of the spectral peek occurs together with reduction in peak power and an increase for off-signal frequencies. A condition that is called spectral leakage.

One notes that the integral value of the spectral energy remains constant at (amp/2)*length(FFT)

print(np.sqrt(np.sum(np.abs(yy0)**2)),np.sqrt(np.sum(np.abs(yy3)**2)))

152.42016636452271 148.52080632743807

yy3=np.fft.rfft(np.real(xx))
xx3=np.fft.irfft(yy3)

# plot signal
fig = plt.figure("figsize",[15,5])
plt.plot(tt,np.real(xx3));
plt.grid(True)
plt.show()

Observation#

The inverse FFT of the complex spectrum generates the original time series, even if the spectrum is spead out over a wider frequency range.

Noise processing#

We now add noise to the sinusoidal signal

plt.rcParams.update({'font.size': 14}) # increase fontsize in figures

fs = 1024   #Hz
nfft=1024
nwin=nfft

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

N = fs
time = np.arange(N) / fs

# signal description
arms=2
amp = arms*np.sqrt(2)
freq = 30

var = 0.001
noise_power = var * fs / 2 
# The factor of fs/2 converts the per sample noise power (variance) 
# to a noise bandwith of -1/T to 1/T (negative and positive frequencies) 

xx = amp*np.sin(2*np.pi*freq*time)
rms= np.sqrt(np.mean(xx**2))

# add noise to signal
np.random.seed(1)
xx += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

print("Signal RMS : ",amp/np.sqrt(2), rms)

#
yy=np.fft.rfft(xx,nfft)
pp=np.abs(yy)**2

ff=np.arange(int(nfft/2)+1)*fs/nfft

ppmn=np.mean(pp[ff>50])
ppmx = np.max(pp)
print("Noise ",ppmn, "Peak ",ppmx)

fig = plt.figure("figure.figsize",[15,10])

plt.subplot(211)
plt.plot(time,xx)
plt.plot(time,amp/np.sqrt(2)+0*time,'r')

plt.subplot(212)
plt.semilogy(ff, pp)
plt.plot(ff,ppmn+0*ff,'r')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Spectrum [V**2]')
plt.grid(True)

Signal RMS :  2.0 2.0
Noise  513.3273885838644 Peak  2073698.6295588852

../_images/9fcd3b9c5c6bf77f49da8dc8c8eee9af191e09eb84363334480f78cb3281498b.png

print(var, 2*ppmn/(N**2))
print(arms, ppmx/(N**2))

0.001 0.0009790942928006446
2 1.9776331229771473

Observation#

To obtain noise spectrum (variance per sample), the mean value of the rfft power spectrum must be scaled by \(2/N^2\) To obtain the RMS of the signal, the peak value of the rfft power spectrum must be scaled by \(1/N^2\)

The extimates are not exactly equal as the noise is random and consequently also all extimates of noisy data are randomvariables.

Spectrum implementation#

For future spectral analysis we can define the following support function that scales the power spectrum accordingly

# initial custom spectrum estimation
def mSpectrum0(xx,nfft,fs,mode=None):
    # xx: real valued time series if at least nfft samples
    # nfft: length of fft in samples
    # fs: sampling frequency
    # mode: RMS for scaling of spectral peakes to signal RMS
    #       VAR for scaling of spectral mean values to per sample variance
    #       None for power spectrum (FFT output) (Default)
    # output: P power estimates according to mode
    #         ff frequency vector
    yy=np.fft.rfft(xx[:nfft],nfft)
    ff=np.arange(int(nfft/2)+1)*fs/nfft

    pp=np.abs(yy)**2
    if mode=="RMS":
        return pp/(nfft**2), ff
    elif mode=="VAR":
        pp[1:] *= 2
        return pp/(nfft**2), ff
    else:
        return pp, ff

P1,ff = mSpectrum0(xx,nfft,fs)
P2,ff = mSpectrum0(xx,nfft,fs,"RMS")
P3,ff = mSpectrum0(xx,nfft,fs,"VAR")


fig = plt.figure("figure.figsize",[15,10])
plt.semilogy(ff,P1);
plt.semilogy(ff,P2);
plt.semilogy(ff,P3);
plt.xlim(0,100);

../_images/f14b3cdf6c40b1fa57e67d25bcab073c9bfdd6db220ec2bd55684af3faee9095.png

Windowing#

The implementation of the FFT assumes that signal is periodic. Most interesting signals are aperiodic, resulting in what is called spectral leakage. The above signal with a frequency of 5.5 is such a case.

Standard remedy is to taper the signal such that both ends of the time series are equal and very often close to zero. A typical window is the Hann window,

(214)#\[\begin{equation} hann(n)=\frac{1}{2}\big(1-\cos(2\pi \frac{n}{N-1})\big) \end{equation}\]

In the following we assume that the window length is equal to FFT length, that is we assume first NO zero padding and analyse later the case of zero padding

win=np.hanning(nwin)

plt.plot(win);
plt.show()
plt.plot(xx*win);
plt.show()

../_images/c000562442dc65b28a3aee52a34044053a45bceb67b915e0ed7d485bf5cea0bf.png

../_images/9d8cbe0f031b11fad1380c7d8c19df8da82008a850ceceb228b3ed52ee8cb95f.png

Observations#

Windowing reduces

the influences of initial and final data
the effective length of the data

Conclusions#

one should expect that

the RMS value of the signal is modified
the noise variance is modified

P11,ff = mSpectrum0(xx*win,nfft,fs)
P12,ff = mSpectrum0(xx*win,nfft,fs,"VAR") # for noise
P13,ff = mSpectrum0(xx*win,nfft,fs,"RMS") # for signal


fig = plt.figure("figure.figsize",[15,10])
plt.semilogy(ff,P11);
plt.semilogy(ff,P12);
plt.semilogy(ff,P13);
plt.xlim(0,100)

print(var, np.mean(P12[50:]))
print(arms, np.max(P13))

0.001 0.00039303184120769764
2 0.5095595110138033

../_images/0a94eec39819bb7a435784fcf6ea4775ac780e0b1b81d23fcfdc1635350ed9a4.png

Observation#

Both, the RMS estimate (peak value) and the per sample Var value (noise) are not correct. The obtained values are too low due to the windowing influence. This can be corrected by considering the window shape.

print(np.mean(win*win))
print(np.mean(win))

print(var, np.mean(P12[50:])/np.mean(win*win))
print(arms, np.max(P13)/win.mean()**2)

3746337890625
49951171875
001 0.0010491094308157248
2.0422248167831363

# updated custom spectrum estimation
def mSpectrum1(xx,win,fs,mode=None):
    # xx: real valued time series if at least nfft samples
    # nfft: length of fft in samples
    # fs: sampling frequency
    # mode: RMS for scaling of spectral peakes to signal RMS
    #       VAR for scaling of spectral mean values to per sample variance
    #       None for power spectrum (FFT output) (Default)
    # output: P power estimates according to mode
    #         ff frequency vector
    nwin=len(win)
    nfft=nwin
    yy=np.fft.rfft(xx[:nwin]*win,nfft)
    ff=np.arange(int(nfft/2)+1)*fs/nfft

    pp=np.abs(yy)**2
    if mode=="RMS":
        return pp/win.sum()**2, ff
    elif mode=="VAR":
        pp[1:] *=2
        return pp/(nwin*(win*win).sum()), ff
    else:
        return pp, ff

P21,ff = mSpectrum1(xx,win,fs)
P22,ff = mSpectrum1(xx,win,fs,"VAR")
P23,ff = mSpectrum1(xx,win,fs,"RMS")


fig = plt.figure("figure.figsize",[15,10])
plt.semilogy(ff,P21);
plt.semilogy(ff,P22);
plt.semilogy(ff,P23);
plt.xlim(0,100)

# test estimates
print(var, np.mean(P22[50:]))
print(arms, np.max(P23))

0.001 0.0010491094308157246
2 2.0422248167831363

../_images/89bd0ee9184e995caf00d3dccfda8b9720f85b1cac4c32412969b514f13fd2f9.png

Observations#

The general scaling functions of spectras are

RMS estimation: 1/win.sum()**2
VAR estimation: 2/(nwin*(win*win).sum()), whereby the factor 2 is only for nonzero frequencies

For rectangular windows (that is no windows) we have win.sum()=nwin and (win*win).sum()=nwin

Power spectral density (PSD)#

The VAR estimation provides the per sample noise variance expressed in V^2/sample . To obtain the PSD that typically is expressed in V^2/Hz one needs to consider the spectral binwidth in terms of the sampling frequency \(f_s\).

Binwidth in Hz is for \(N_{fft}= N_{win}\) given by

(215)#\[\begin{equation} \Delta f = \frac{f_s}{N_{win}} \end{equation}\]

and therfore the PSD related scaling of the spectra

PSD estimation: 2/(nwin*(win*win).sum())*nwin/fs or 2/(fs*(win*win).sum())

Zero padding#

So far, we assumed that the length of the fft is equal to the length of the window. In the following we allow the fft length to be longer than the window length, wherby the added samples are kept at zero (zero padding).

# case 1
nfft1=nwin
uu=xx[:nwin]*win
yy1=np.fft.rfft(uu,nfft1)
ff1=np.arange(int(nfft1/2)+1)*fs/nfft1

pp1=np.abs(yy1)**2

# case 2
nfft2=4*nwin
uu=np.zeros(nfft2)
uu[:nwin] = xx[:nwin]*win
yy2=np.fft.rfft(uu,nfft2)
ff2=np.arange(int(nfft2/2)+1)*fs/nfft2

pp2=np.abs(yy2)**2
    
fig = plt.figure("figure.figsize",[15,10])
plt.semilogy(ff1,pp1);
plt.semilogy(ff2,pp2);

print(np.max(pp1),np.max(pp2))
print(np.mean(pp1[ff1>40]),np.mean(pp2[ff2>40]))

534311.8738208098 534311.8738208096
205.5119088000577 205.5145016630521

../_images/cbc2b47ae67bf53281e53606ad399d89b3e285c2d1f535ed7468aa9af674952d.png

Observations#

We note that the unscaled spectra estimates do not vary when adding zeros to the data before FFT (zero padding). All what is changing is the number of frequency bins to cover the whole spectrum and consequently, zero padded FFTs are also called interpolationg FFT.

One should expect that scaling for the RMS estimate (peak of spectrogram) does not depend on length of zero padded FFT but only on length of window

Open question: does zero padding change the estimate of the PSD?#

The PSD is given as the time averaged spectral energy

(216)#\[\begin{equation} P(f) = \frac{1}{T}\Big|\int_{0}^{T}x(t)\exp{(-2\pi i f t)} dt\Big|^2 \end{equation}\]

or for discrete sampling

(217)#\[\begin{equation} P(f) = \frac{1}{N\Delta t}\Big|\sum_{n=0}^{N-1}x(n)\exp{(-2\pi i f n\Delta t)} \Delta t\Big|^2 \end{equation}\]

where \(N\) is the number of samples of the signal and \(\Delta t\) is the sampling interval, i.e. \(N\Delta t=T\) is the length of the signal in seconds and \(n\Delta t\) is the time when the n-th samples was measured.

The quantity \(P(f)\) is really a time averaged signal energy

Assume zero padding from \(N\) to \(K>N\), the sum does not change as \(x(n)=0\) for all \(n>N\), but the frequency interpretation does:

Assume \(f_k =k(f_s/K)\) for \(k=0,1, ...,K-1\) with \(f_s = 1/\Delta t\)

(218)#\[\begin{equation} P(k) = P(\frac{kf_s}{K}) = \frac{1}{f_sN}\Big|\sum_{n=0}^{N-1}x(n)\exp{(-2\pi i kn/K)}\Big|^2 \end{equation}\]

The PSD is independent of the DFT size K, but it only depends on the number of time domain samples N

Comment#

Zero padding is only useful if one is interested in the interpolating aspect of the zero-padded FFT, e.g. for RMS estimates of tonal signals. For broadband PSD estimation zero-padded FFT seems of little advantage.

Custom Spectrum implementation#

# custom spectrum estimation
def mSpectrum(xx,win,fs,nfft,mode=None):
    # xx: real valued time series if at least nfft samples
    # nfft: length of fft in samples
    # fs: sampling frequency
    # mode: RMS for scaling of spectral peakes to signal RMS
    #       VAR for scaling of spectral mean values to per sample variance
    #       PSD for scaling of power spectral density to (v^2/Hz)
    #       None for power spectrum (FFT output) (Default)
    # output: P power estimates according to mode
    #         ff frequency vector
    nwin=len(win)
    if nfft<nwin: nfft=nwin # force nfft to be at least window length
    nfft=int(2**np.ceil(np.log2(nfft))) # increase nfft to next power of 2
    #
    # copy data to working buffer zero-pad if required
    uu=np.zeros(nfft)
    uu[:nwin]=xx[:nwin]*win
    
    yy=np.fft.rfft(uu,nfft)
    ff=np.arange(int(nfft/2)+1)*fs/nfft

    pp=np.abs(yy)**2
    if mode=="RMS":
        return pp/win.sum()**2, ff
    elif mode=="VAR":
        pp[1:] *=2
        return pp/(nwin*(win*win).sum()), ff
    elif mode=="PSD":
        pp[1:] *=2
        return pp/(fs*(win*win).sum()), ff
    else:
        return pp, ff

# case 1
nfft1=nwin
# case 2
nfft2=4*nwin

#----------------------------------------------
pp1,ff1=mSpectrum(xx,win,fs,nfft1,"RMS")
pp2,ff2=mSpectrum(xx,win,fs,nfft2,"RMS")

fig = plt.figure("figure.figsize",[15,10])
plt.semilogy(ff1,pp1);
plt.semilogy(ff2,pp2);
plt.show()

print("RMS: ",np.max(pp1),np.max(pp2))

#----------------------------------------------
pp1,ff1=mSpectrum(xx,win,fs,nfft1,"PSD")
pp2,ff2=mSpectrum(xx,win,fs,nfft2,"PSD")

fig = plt.figure("figure.figsize",[15,10])
plt.semilogy(ff1,pp1);
plt.semilogy(ff2,pp2);
plt.show()

print("PSD: ", np.mean(pp1[ff1>40]),np.mean(pp2[ff2>40]))

../_images/360f4612603a8d93098c8f99d8523ea59abfec5cca1847497f212b6d2a89ec43.png

RMS:  2.0422248167831363 2.042224816783136

../_images/2263f7be5978a25c92b74e53c57c7f252d95e544659e560245b073898bb2a050.png

PSD:  0.0010463094086024442 0.0010463226094770898

BLM_1: Basic signal spectral analysis

Contents

BLM_1: Basic signal spectral analysis#

Observations#

Complex numbers#

Trigonometry#

Pythagoras#

Problem#

Hint#

Solution#

Definitions#

complex number (real + i*imaginary)#

conjugate of a complex number#

absolute squared complex value by signal multiplied with its conjugate#

Alternative method#

General frequencies#

Frequency analysis#

Discrete Fourier transform (DFT)#

Observation#

Fast Fourier Transform#

Comments#

FFT implementation#

Observation#

General signal frequency#

Observations#

Conclusions#

Spectral leakeage#

Observation#

Noise processing#

Observation#

Spectrum implementation#

Windowing#

Observations#

Conclusions#

Observation#

Observations#

Power spectral density (PSD)#

Zero padding#

Observations#

Open question: does zero padding change the estimate of the PSD?#

Comment#

Custom Spectrum implementation#