Spectral Windowing and Variance Reduction

This script illustrates how to use windowing to get localized signal information and how to reduce the spatial/temporal variance. We will be using some of the simplifications outlined in the intro_spectral1d.m script.

Make a grid and wavenumbers
Sample Function
PSD with no windowing
Windowing
Variance Reduction

close all;clear;

Make a grid and wavenumbers

N = 1024;
L = 20;
z = linspace(-1,1,N+1);
z = z(1:end-1);
x = z*L;
dx = x(2)-x(1);
dk = pi/L;
ksvec = 0:N/2-1;
k = ksvec'*dk;

Sample Function

For this example we will use the a more complicated function, the Jacobi Elliptic function. This function has clearly a clearly defined spectrum with peaks at 1.885 and 5.498. We then add a random perturbation to this signal.

[fs,~,~]=ellipj(2*pi*x/2.5,0.75);
f1=fs+0.005*randn(size(x));
fft_f1 = abs(fft(f1));

PSD with no windowing

Without any windowing we can see that the spectrum is extremely problematic. While there are peaks at the correct wavenumbers, there is power at every wavenumber and there are small oscillations at large wavenumbers.

figure(1)
subplot(2,1,1)
plot(x,f1)
xlabel('x')
title('Sample function and its raw spectrum')
subplot(2,1,2)
semilogy(k,fft_f1(1:end/2))
xlabel('k')
ylabel('Raw PSD')

Windowing

These problems can be partially solved by windowing. Windowing is multiplying the original function by a "window" function that is localized You can think of windowing like making a packet from your signal. Here we specify the center point (N/2) and then the windowing function linearly decreases from there.

No2 = N/2;
mywin = [0:No2-1 No2:-1:1]/No2;
f1w = mywin.*f1;
fft_f1w = abs(fft(f1w));
figure(2)
subplot(2,1,1)
plot(x,f1,'b',x,f1w,'r')
title('Sample function, windowed function and their spectra')
subplot(2,1,2)
semilogy(k,fft_f1(1:end/2),'b',k,fft_f1w(1:end/2),'r')

You can see that this has decreased the power at most wavenumbers however we still have problems. The windowed spectrum is more choppy.

Variance Reduction

Using abs(fft(f)) gives an estimate of the spectrum with a lot of variance (technically an astonishing 100%). In practice things are rarely that bad, but when a signal is long the variance can be reduced by splitting the signal up, windowing each piece and summing.

No4 = N/4;
mywin = [0:No4/2-1 No4/2:-1:1]/(No4/2);
myspec = zeros(No4,1)';
for ii=0:3
    fnow = mywin.*f1(ii*No4+1:(ii+1)*No4);
    myspec = myspec+abs(fft(fnow));
end
dk4 = 4*pi/L;
k4svec = 0:No4/2-1;
k4 = k4svec'*dk4;
figure(3)
semilogy(k4,myspec(1:end/2),'k',k,fft_f1w(1:end/2),'r')
title('Windowed and variance reduced spectra')
legend('Smoothed','Original')
xlabel('k')
ylabel('PSD')
axis([0 20 2e-3 5e2])
grid on

The other way of reducing variance is by just directly smoothing the spectrum.

spec0sm = fft_f1w;
for ii=3:N-2
    spec0sm(ii) = sum(fft_f1w(ii-2:ii+2))/5;
end
figure(4)
semilogy(k,spec0sm(1:end/2),'g',k,fft_f1w(1:end/2),'r')
title('Windowed and variance reduced spectra')
legend('Smoothed','Original')
axis([0 20 2e-3 5e2])