LPFCSD_Example_3

Ivan Selesnick,
Polytechnic Institute of NYU,
August 2013

Contents

Start

clear
close all

printme = @(filename) print('-dpdf', sprintf('figures/LPFCSD_Example_3_%s', filename));

set(0, 'DefaultAxesFontSize', 8)

Load data

y = load('wl1_col02.txt');

Fs = 6.25;          % Fs : sampling rate (samples/second)
N = length(y);
n = 1:N;
t = n/Fs;

figure(1)
clf
plot(t, y);
xlim([0 N/Fs])
ax1 = axis;
axis(ax1);
title('Data')

Run LPF/CSD algorithm

% Set filter parameters
fc = 0.05;       % fc : cut-off frequency (cycles/sample)
d = 1;           % d : filter order parameter (d = 1, 2, or 3)

lam0 = 0.005;
lam1 = 1;

Nit = 50;
% Nit : number of iterations

mu = 0.2;
% mu : ADMM parameter for LPF/CSD

% Run LPF/CSD algorithm
[x, f, cost] = lpfcsd(y, d, fc, lam0, lam1, Nit, mu);

txt = sprintf('fc = %.2e, lam0 = %.2e, lam1 = %.2e', fc, lam0, lam1);

Display cost function history

figure(4)
clf
plot(cost)
title('Cost function history')
xlabel('Iteration')

Run-time calculation

if 0
    % Run-time computation
    tic
    for ii = 1:50
        [x_, f_, cost_] = lpfcsd(y, d, fc, lam0, lam1, Nit, mu);
    end
    time_lpfcsd = toc/50;
    fprintf('LPF/CSD run-time: %.4f seconds (length = %d samples, %d iterations)\n', time_lpfcsd, length(y), Nit);
end
figure(3)
clf
subplot(3, 1, 1)
plot(t, y)
xlim([0 N/Fs])
ax = axis;
title('Data')
axis(ax)

subplot(3, 1, 2)
plot(t, x)
xlim([0 N/Fs])
title('CSD component [Output of LPF/CSD algorithm]')

subplot(3, 1, 3)
plot(t, y - x)
xlim([0 N/Fs])
xlabel(txt)
title('Data with CSD component subtracted')
xlabel('Time (sec)')

orient landscape
printme('lpfcsd')

Detail

xlim1 = 180+[0 100];

subplot(3, 1, 1)
title('Data (detail)')
xlim(xlim1)

subplot(3, 1, 2)
title('CSD component [Output of LPF/CSD algorithm] (detail)')
xlim(xlim1)

subplot(3, 1, 3)
title('Data with CSD component subtracted (detail)')
xlim(xlim1)

printme('lpfcsd_detail')

Periodograms

Z = pwelch(y - x);

L = length(Z)-1;
f = (0:L)/L * Fs/2;

Y = pwelch(y);

figure(2)
clf
plot(f, 20*log10(Y), f, 20*log10(Z),'r');
xlabel('Frequency (Hz)')
title('Welch periodogram')
legend('Data', 'Data with LPF/CSD')
xlim([0 Fs/2])

printme('pwelch')

Print figure

figure(102)
clf

subplot(5, 1, 1)
plot(t, y, 'blue')
xlim([0 N/Fs])
ax = axis;
title('Original data');
axis(ax)
xlabel('Time (sec)')
ylim([-25 75])

subplot(5, 1, 2)
plot(t, x, 'black')
xlim([0 N/Fs])
title('CSD component (x)');
xlabel('Time (sec)')
ylim([-25 75])


subplot(5, 1, 3)
plot(t, y - x, 'red')
xlim([0 N/Fs])
title('Corrected data (data with CSD component subtracted)');
xlabel('Time (sec)')

subplot(5, 1, [4 5])
plot(f, 20*log10(Y), 'color', 'blue');
line(f, 20*log10(Z), 'color', 'red');
legend('Original data', 'Corrected data', 'location', 'northeast')
xlabel('Frequency (Hz)')
title('Welch periodogram (dB)')
xlim([0 Fs/2])
ylim([-40 50])

orient tall
printme('output')
set(0, 'DefaultAxesFontSize', 'remove')