Using iterative hard thresholding to recreate a signal. dtw and idwt can easily be replaced with fft and ifft.

Iterative Hard Thresholding

    clear all
    close all
    clc
    % Load in an image
    I = double(imread('~/Desktop/not-used-frequently/pictures_for_project/lenna.jpg'));
    I2 = mean(I,3);
    I3 = imresize(I2,[512,512]);

    % Take Measurements (pixels)
    sz = size(I3);
    n = sz(1)*sz(2);
    p = .3;
    rp = randperm(n);

    % This line is the image --> pixel samples operation
    % implementing the linear function A(.)
    y = I3(rp(1:floor(p*n)));

    % This is the operation that takes a vector of samples
    % and forms the full image, placing the samples at the right places
    % and zeros elsewhere. This is the adjoint of A(.)
    ys = zeros(size(I3));
    ys(rp(1:floor(p*n))) = y;


    % Display sampled image
    subplot(1,2,1)
    imagesc(I3); colormap gray
    subplot(1,2,2)
    imagesc(ys); colormap gray
    drawnow

    % Now let's implement IHT

    % Note that we're going to go after the wavelet transform coefficients...
    % So, my effective observation model is

    % y = A(T^(-1)(T(z)))
    % let T(z) = x, then
    % y = A(T^(-1)(x))
    % let phi = A(T^(-1)(.))
    % our model is y = A(x)
    iterations = 200;
    s = 5000;
    %h = daubcqf(2);
    % initialize
    xold = zeros(size(I3));
    for i=1:iterations
        % First compute phi(xold)
        t1 = idwt2_full(xold);
        temp = t1(rp(1:floor(p*n)));

        % Now compute y-phi(xold)
        temp2 = y-temp;

        % Now compute phi*((y-phi(xold))
        % (first "blow this up" into a full matrix)
        temp3 = zeros(size(I3));
        temp3(rp(1:floor(p*n))) = temp2;
        % (then take the forward wavelet transform)
        temp3 = dwt2_full(temp3);

        % Now compute xold + phi*((y-phi(xold))
        temp4 = xold + temp3;

        % Finally do hard thresholding
        s = round(5000 + i*n/1000);
        if s > n s = n; end
        s = 5000;
        
        tt = reshape(temp4,n,1);
        [srt,inx] = sort(abs(tt),'descend');
        temp5 = zeros(size(temp4));
        temp5(inx(1:s)) = tt(inx(1:s));
        abs(tt(inx(s)))
        
        %ma = max(tt);
        %temp5 = zeros(sz);
        %cut_off = 0.01 * ma;
        %i = tt < ma;
        %temp5(i) = 0;

        % At each step, plot the current estimate...
        imagesc(idwt2_full(temp5)); colormap gray
        title(num2str(i))
        drawnow

        % update
        xold = temp5;
    end
    imagesc(idwt2_full(temp5)); colormap gray
        title(num2str(i))
        drawnow