Welch's Method with Windows

As usual, the purpose of the window function $ w(n)$ (Chapter 3) is to reduce side-lobe level in the spectral density estimate, at the expense of frequency resolution, exactly as in the case of sinusoidal spectrum analysis.

When using a non-rectangular analysis window, the window hop-size $ R$ cannot exceed half the frame length $ M/2$ without introducing a non-uniform sensitivity to the data $ x(n)$ over time. In the rectangular window case, we can set $ R=M$ and have non-overlapping windows. For Hamming, Hanning, and any other generalized Hamming window, one normally sees $ R=(M-1)/2$ for odd-length windows. For the Blackman window, $ R\approx M/3$ is typical. In general, the hop size $ R$ is chosen so that the analysis window $ w$ overlaps and adds to a constant at that hop size. This consideration is explored more fully in Chapter 8. An equivalent parameter used by most matlab implementations is the overlap parameter $ M-R$ .

Note that the number of blocks averaged in (6.27) increases as $ R$ decreases. If $ N_x\geq M$ denotes the total number of signal samples available, then the number of complete blocks available for averaging may be computed as

$\displaystyle K = \left\lfloor \frac{N_x-M}{R}\right\rfloor + 1,$ (7.28)

where the floor function $ \left\lfloor x\right\rfloor $ denotes the largest integer less than or equal to $ x$ .


Matlab for Welch's Method

Octave and the Matlab Signal Processing Toolbox have a pwelch function. Note that these functions also provide confidence intervals (not covered here). Matlab for generating the data shown in Fig.6.1 (employing Welch's method) appears in Fig.6.2.

Figure: Matlab listing for generating data plotted in Fig.6.1.

 
M = 32;
Ks=[1 8 32 128]
nkases = length(Ks);
for kase = 1:4
  K = Ks(kase);
  N = M*K;
  Nfft = 2*M; % zero pad for acyclic autocorrelation
  Sv = zeros(Nfft,1); % PSD "accumulator"
  for m=1:K
    v = randn(M,1);  % noise sample
    V = fft(v,Nfft);
    Vms = abs(V).^2; % same as conj(V) .* V
    Sv = Sv + Vms;   % sum scaled periodograms
  end
  Sv = Sv/K;     % average of all scaled periodograms
  rv = ifft(Sv); % Average Bartlett-windowed sample autocor.
  rvup = [rv(Nfft-M+1:Nfft)',rv(1:M)']; % unpack FFT
  rvup = rvup/M; % Normalize for no bias at lag 0
end


Next Section:
Filtered White Noise
Previous Section:
Welch's Method