Example: Basis pursuit denoising (BPD) with the STFT
Estimation of pulse in white Gaussian noise. The STFT is used.
Contents
Misc
clear
close all
I = sqrt(-1);
dB = @(x) 20 * log10(abs(x));
Make signal
N = 100; % N : signal length n = 0:N-1; f = 0.15; K = 20; % K : duration of pulse M = 55; % M : start-time of pulse A = 3; s = zeros(size(n)); s(M+(1:K)) = A*sin(2*pi*f*(1:K)) .* hamming(K)'; % s : pulse figure(1) clf plot(n, s) title('Noise-free signal [s]')
Define transform
R = 16; % R : frame length Nfft = 32; % Nfft : FFT length in STFT (Nfft >= R) [A, AH, normA] = MakeTransforms('STFT', N, [R Nfft]);
Reconstruction error = 0.000000 Energy ratio = 1.000000
Validate normA calculation
% imp : an impulse in transform domain imp = AH(zeros(1, N)); n0 = 10; k0 = 5; imp(n0, k0) = 1; [normA norm(A(imp))] % should be same
ans = 0.5000 0.5000
Make noisy data
randn('state', 0); % set state for reproducibility of example sigma = 0.5; y = s + sigma * randn(size(s));
Display noisy data and its STFT
figure(1) clf subplot(2, 1, 1) plot(n, y) title('Noisy data [y]') AHy = AH(y); tt = R/2 * ( (0:size(AHy, 2)) - 1 ); % tt : time axis for STFT subplot(2, 1, 2) imagesc(tt, [0 0.5], dB(AHy(1:Nfft/2+1, :)), max(dB(AHy(:))) + [-50 0]) axis xy xlim([0 N]) ylim([0 0.5]) title('STFT of Noisy data') ylabel('Frequency') xlabel('Time') cm = colormap(gray); % cm : colormap for STFT cm = cm(end:-1:1, :); colormap(cm)
Solve BPD
beta = 2.5; lam = beta * sigma * normA; mu = 0.1; Nit = 100; [c, cost] = BPD(y, A, AH, lam, mu, Nit); x = A(c); % Display cost function history to observe convergence of algorithm. figure(1) clf plot(cost) title('Cost function history') xlabel('Iteration')
Display denoised signal
cmax = max(abs(c(:))); figure(1) clf subplot(2, 1, 1) plot(n, real(x)) title('Denoised signal') ylabel('real( Ac )') ch = colorbar; set(ch, 'visible', 'off') subplot(2, 1, 2) imagesc(tt, [0 0.5], dB(c(1:Nfft/2+1, :)), dB(cmax) + [-50 0]) % imagesc(tt, [0 1], dB(c), dB(cmax) + [-50 0]) axis xy xlim([0 N]) ylim([0 0.5]) % title('BPD STFT coefficients'); title(sprintf('BPD STFT coefficients (beta = %.2f)', beta)); ylabel('Frequency') xlabel('Time') colormap(cm) ch = colorbar;
Optimality scatter plot
I = sqrt(-1); g = (1/lam) * AH(y - A(c)) .* exp(-I*angle(c)); figure(2) clf plot3( real(g(:)), imag(g(:)), abs(c(:)), '.') zlabel('abs( c )') xlabel('Re(A^H(y - A c) e^{-j\anglec})/\lambda') ylabel('Im(A^H(y - A c) e^{-j\anglec})/\lambda') title('Optimality scatter plot') grid off theta = linspace(0, 2*pi, 200); line( cos(theta), sin(theta), 'color', [1 1 1]*0.5) line([1 1], [0 0], [0 cmax], 'color', [1 1 1]*0.5) xm = 1.2; xlim([-1 1]*xm) ylim([-1 1]*xm) line([-1 1]*xm, [1 1]*xm, 'color', 'black') line([1 1]*xm, [-1 1]*xm, 'color', 'black') % Animate: vary view orientation az = -36; for el = [40:-1:5] view([az el]) drawnow end for az = [-36:-3] view([az el]) drawnow end
Save figure to file
MyGraphPrefs('on') figure(1) clf subplot(4, 1, 1) line(n, y) mytitle('Noisy data'); ylabel('y') xlabel('Time') ylim([-3 3]) subplot(4, 1, 2) imagesc(tt, [0 0.5], dB(AHy(1:Nfft/2+1, :)), max(dB(AHy(:))) + [-50 0]) % imagesc(tt, [0 1], dB(abs(AH(y))), dB(cmax) + [-50 0]) axis xy colormap(cm) xlim([0 N]) ylim([0 0.5]) mytitle('STFT of noisy data'); xlabel('Time') ylabel('Frequency') box off subplot(4, 1, 3) imagesc(tt, [0 0.5], dB(c(1:Nfft/2+1, :)), dB(cmax) + [-50 0]) % imagesc(tt, [0 1], dB(c), dB(cmax) + [-50 0]) axis xy colormap(cm) xlim([0 N]) ylim([0 0.5]) mytitle(sprintf('BPD STFT coefficients (beta = %.2f)', beta)); xlabel('Time') ylabel('Frequency') box off subplot(4, 1, 4) line(n, real(A(c))) mytitle('Denoised signal'); ylabel('real( Ac )') xlabel('Time') ylim([-3 3]) % print figure to pdf file orient tall print -dpdf figures/BPD_example_stft % print figure to eps file set(gcf, 'PaperPosition', [1 1 4 7.5]) print -deps figures_eps/BPD_example_stft MyGraphPrefs('off')
Animate: vary lambda
for beta = 0.1:0.1:3 lam = beta * sigma * normA; [c, cost] = BPD(y, A, AH, lam, mu, Nit); x = A(c); figure(10) clf subplot(2, 1, 1) imagesc(tt, [0 0.5], dB(c(1:Nfft/2+1, :)), dB(cmax) + [-50 0]) % imagesc(tt, [0 1], dB(abs(c)), [-20 20]) axis xy xlim([0 N]) ylim([0 0.5]) title(sprintf('BPD STFT coefficients : \\lambda = \\beta\\sigma ||a||, \\beta = %.2f', beta)); ylabel('Frequency') xlabel('Time') colormap(cm) subplot(2, 1, 2) plot(n, real(x)) ylim([-1 1]*4) title(sprintf('Denoised signal, \\beta = %.2f', beta)); ylabel('real( Ac )') xlabel('Time') drawnow end