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')