demo_FFT_2
Complex data example
Ivan Selesnick March 2021
Contents
Start
clear set_plot_defaults % Set random number generator for reproducibility rng('default') rng(1) rmse = @(err) sqrt(mean(abs(err(:)).^2));
Define transform
The transform A is wide-matrix inverse DFT. A is of size 100 x 256. A satisfies A A' = I
M = 100; % M: signal length N = 256; % N: transform length truncate = @(x, M) x(1:M); % A: inverse DFT (normalized for Parsevals property) AH = @(x) fft(x, N)/sqrt(N); A = @(X) truncate(ifft(X), M) * sqrt(N); normA = sqrt(M/N);
Verify transform properties
Verify properties of A
% Verify that A AH = I y = randn(M, 1); err = A(AH(y)) - y; max(abs(err)) % Verify A as a matrix Amtx = dftmtx(N)'/sqrt(N); Amtx = Amtx(1:M, :); x = randn(M, 1); err = Amtx' * x - AH(x); max(abs(err))
ans = 4.4409e-16 ans = 4.5776e-16
Create data
Sinusoids plus white noise
f1 = 0.1; f2 = 0.22; m = (0:M-1)'; % s : true signal C1 = 2*exp(1j * pi/3); C2 = 1*exp(-1j * pi/4); s = C1 * exp(1j * 2*pi * f1*m) + C2 * exp(1j * 2*pi* f2*m); sigma = 1; % sigma : noise standard deviation w = randn(M, 1) + 1j * randn(M, 1); % w : white Gaussian noise y = s + sigma * w; % y : signal plus noise figure(1) clf subplot(2, 1, 1) plot(m, real(s), m, imag(s)) title('Noise-free data') xlabel('Time') subplot(2, 1, 2) plot(m, real(y), m, imag(y)) title('Noisy data') xlabel('Time')

L1 norm
Least squares with L1 norm regularization
lam_L1 = 1.0; c_L1 = sparseLS_L1(y, A, AH, 1, lam_L1); x_L1 = A(c_L1); rmse_L1 = rmse(x_L1 - s) figure(1) clf % Signals subplot(2, 1, 1) plot(m, real(s), m, real(x_L1)) title(sprintf('Denoising using L1 norm, RMSE = %.3f', rmse(x_L1 - s))) legend('True signal', 'L1 norm'); ylim1 = [-1 1]*max(abs(y)) * 1.1; ylim(ylim1) xlabel('Time (n)') % Sparse coefficients subplot(2, 1, 2) plot(1:N, abs(c_L1), 'o') title('Sparse Fourier coefficients (L1 norm)') xlim([0 N]) xlabel('Frequency (index)') ylim([0 15]) drawnow print -dpdf -bestfit figures/demo_FFT_2_L1norm
rmse_L1 = 0.5297

GMC penalty
Least squares with GMC regularization
lam_GMC = 2.0; gamma = 0.8; c_GMC_ = sparseLS_GMC(y, A, AH, 1, lam_GMC, gamma);
Altnerate algorithm
c_GMC = sparseLS_GMC_ver2(y, A, AH, 1, lam_GMC, gamma); figure(2) clf % plot(1:N, abs(c_GMC), 1:N, abs(c_GMC_ver2)) plot( 1:N, abs(c_GMC_), 1:N, abs(c_GMC_ - c_GMC), 'o') legend('coefficients', 'erorr') x_GMC = A(c_GMC); rmse_GMC = rmse(x_GMC - s)
rmse_GMC = 0.3395

figure(1) clf % Signals subplot(2, 1, 1) plot(m, real(s), m, real(x_GMC)) title(sprintf('Denoising using GMC, RMSE = %.3f', rmse(x_GMC - s))) legend('True signal', 'GMC penalty'); ylim(ylim1) xlabel('Time (n)') % Sparse coefficients subplot(2, 1, 2) plot(1:N, abs(c_GMC), 'o') title('Sparse Fourier coefficients (GMC penalty)') xlim([0 N]) xlabel('Frequency (index)') ylim([0 15]) drawnow print -dpdf -bestfit figures/demo_FFT_2_GMC

Plot sparse coefficients
figure(1) clf plot(1:N, abs(c_GMC), 'o', 1:N, abs(c_L1), 'x') legend('GMC', 'L1', 'location', 'north') title('Sparse Fourier coefficients') xlim([0 N])

Plot solutions
f = (0:N-1)/N; figure(2) clf % Noise-free signal subplot(2, 2, 1) plot(m, real(s), m, imag(s)) title('Noise-free signal'); ylim1 = [-1 1]*max(abs(y)) * 1.1; ylim(ylim1) xlabel('Time (n)') % Noisy data subplot(2, 2, 2) plot(m, real(y), m, imag(y)) title(sprintf('Noisy signal, \\sigma = %.2f', sigma)); ylim(ylim1) xlabel('Time (n)') subplot(2, 2, 3) plot(f, abs(AH(y))) title('FFT of noisy signal'); xlabel('Frequency') % xlim([0 0.5]) ylim([0 8]) % Sparse coefficients subplot(2, 2, 4) k_L1 = abs(c_L1) > eps; k_GMC = abs(c_GMC) > eps; ph_L1 = plot(f(k_L1), abs(c_L1(k_L1)), 'x'); hold on ph_GMC = plot(f(k_GMC), abs(c_GMC(k_GMC)), 'o', 'markersize', 4); legend([ph_GMC, ph_L1], 'GMC', 'L1 norm') title('Optimized coefficients'); xlabel('Frequency') ylim([0 15]) % xlim([0 0.5]) drawnow orient landscape print -dpdf -fillpage figures/demo_FFT_2
