Building whitening filters

Whitening filters can be constructed from power spectral density (PSD) of noise only data series. The procedure can be summarized in few steps:

  1. Get noise only data
  2. Etimate PSD
  3. Obtain a smoothed version of your PSD estimation
  4. Use the smoothed PSD to calculate the whitening filter

The process starts with the doenload of noise only data.

    noise = matrix('noise.mat');

Assuming the noise is stationary we calculate whitening filter from the noise only part of the first experiment.

    [n1, n2] = noise(1).unpack;

Whitening filter will be calculated from the noise spectra. We assume that the correlations between the two data series are so low that they cannot be properly identified by a cpsd estimation. Therefore we calculate two independent whitening filters for the two output channels. This is equivalent to assume uncorrelated data.

    n1x = n1.lpsd;
    n2x = n2.lpsd;
We are going to split in order to remove low frequency window corrupted data-points,
    n1xx = split(n1x, plist('samples', [6 numel(n1x.y)]));
    n2xx = split(n2x, plist('samples', [6 numel(n2x.y)]));
and finally plot to check results.
    iplot(n1xx, n2xx)

Since the amount of data is limited, we do not properly cover low frequency with the estimated spectrum. We need to exted low frequency part in order to avoid that the fit would introduce unwanted features in a sensible frequency range (e.g. if we skip this part we end-up with a notch on channel 1 that will be clearly visible in the whitened time series).

    xt = [logspace(-5, log10(4.5e-4), 60)]';
    level1 = n1xx.y(1) .* 0.90;
    yt1 = ones(numel(xt),1) .* level1;

    level2 = n2xx.y(1) .* 3.00;
    yt2 = ones(numel(xt),1) .* level2;
and generate an AO from them.
    yunit  = n1xx.yunits;
    fs     = n1xx.fs;    
    ext1   = ao(plist('XVALS', xt, 'YVALS', yt1, 'TYPE', 'FSDATA', 'YUNITS', yunit, 'FS', fs));
    ext2   = ao(plist('XVALS', xt, 'YVALS', yt2, 'TYPE', 'FSDATA', 'YUNITS', yunit, 'FS', fs));
Original data and extension can be joint now.
    nn1xx = join(ext1, n1xx);
    nn2xx = join(ext2, n2xx);
The resulting spectrum che be checked with a plot.
    iplot(nn1xx, nn2xx)

In order to obtain a smooth version of the calculated PSD we perform a frequency domain fit with the 'ao' class method 'sDomainFit', which output a 'parfrac' object that can be used to define the input for our whitener.

      plfit = plist(...
                'AUTOSEARCH', 'on', ...
                'STARTPOLESOPT', 'clog', ...
                'maxiter', 30, ...
                'minorder', 4, ...
                'fittol', 1e-2, ...
                'msevartol', 5e-1, ...
                'plot', 'off');
      s1 = sDomainFit(nn1xx,plfit);
      plfit = plist(...
                'AUTOSEARCH', 'on', ...
                'STARTPOLESOPT', 'clog', ...
                'maxiter', 30, ...
                'minorder', 3, ...
                'fittol', 1e-2, ...
                'msevartol', 5e-1, ...
                'plot', 'off');
      s2 = sDomainFit(nn2xx,plfit);

Frequency series 'aos' can be constructed from the 'parfrac' objects. They will be used as input to 'buildWhitener'.

    plr  = plist('f',logspace(-5, log10(5), 5000)); % define plist with frequency grid

    rs1  = s1.resp(plr); % get frequency response from parfarc object
    mod1 = abs(rs1); % this is a model for a psd, so it must be real

    rs2  = s2.resp(plr); % get frequency response from parfarc object
    mod2 = abs(rs2); % this is a model for a psd, so it must be real


We are now ready to build the whitening filters for the first and differential channels. 'buildWhitener', build a whitening filter from a frequency series input and a 'plist'. The frequency series 'ao' should be representative of the one-sided psd of the noise we want to whiten. 'buildWhitener' perform a fit in the z-domain to the data in order to identify the coefficients (poles and residues) of a digital filter.

    plw = plist(...
          'fs', fs,...
          'FITTOL', 1e-3);

    w1 = buildWhitener1D(mod1, plw);
    w2 = buildWhitener1D(mod2, plw);
If you want more information about 'buildWhitener options, then type
    help ao/buildWhitener1D
in MATLAB command window.

Once we have the whitening filters, we can check their frequency response. It is good practice to ceck filters response on a frequency grid tighter than that used for buildind the filter.

    plr = plist('f', logspace(-5, log10(5), 10000));

    rw1 = w1.resp(plr);
    rw2 = w2.resp(plr);

It is also interesting to check the behavior of the whitening filters on the noise time series.
    nw1 = filter(n1, w1);
    nw2 = filter(n2, w2);

    iplot(nw1, nw2)
    nw1xx = lpsd(nw1);
    nw2xx = lpsd(nw2);

    iplot(nw1xx, nw2xx)
As can be seen, the whitening is not perfect expecially at low frequency but enough to perform a meaningful fit.

Once we are happy with our whitening filters, we can pack them in a matrix (diagonal) and save it. In order to have a diagonal filter matrix we will put empty filter objects out of the diagonal.

    % out of diagonal filter
    fil = filterbank(miir());

    % build a matrix with whitening filters elements
    wf = matrix(w1, fil, fil, w2, plist('shape', [2 2]));'whitening_filter.mat');

©LTP Team