www.gusucode.com > 压缩感知重构算法 压缩感知源码程序 > StOMP_duibi.m

    clc;
clear all;
close all;
sig01=load('cslunwen3.txt');
sig0=sig01;
xfs=50;
points=2048;
f = xfs*(-points:points-1)/points;
figure;
plot(f,sig0);
grid on;
M = 2^11;
qmf = MakeONFilter('Haar',1);
SNR =30;
y=awgn(sig0,SNR);

figure(1);
y1=ifft(y,2*M);
plot(0:length(y1)-1,abs(y1));
grid on;
f = xfs*(-points:points-1)/points;

figure(2);
y2=fft(y1,2*M);
plot(f,y2);
axis([-50 50 -42 0]);
% axis([0 20 -0.1 1.1]);
xlabel('f[MHz]');
ylabel('Normalized PSD');
title('Nyquist rate PSD');
grid on;
for N=48:100:M
A = MatrixEnsemble(N,M);
AM=MatrixEnsemble(2*N,2*M);
L =3;
alpha0 = FWT_PO(y, L, qmf);
Phi=AM;
% generate the vector S
S = Phi * alpha0;
% Solve the StOMP problem
alpha = SolveStOMP(Phi,S,2*M);
sig = IWT_PO(alpha,L,qmf);
mse(N)=twonorm(sig-sig0)/twonorm(sig0);
if N==948
    figure(3);
    subplot(3,1,1);
    plot(f,y); 
    axis([-50 50 -42 0]);
    xlabel('f[MHz]');
    ylabel('PSD');
    title('Nyquist rate PSD');
    % title('(a) Noisy Signal,SNR=15');
    grid on;
    subplot(3,1,2);
    plot(f,sig); 
    axis([-50 50 -42 0]);
    % axis([0 M -1 7]); 
    xlabel('f[MHz]');
    ylabel('PSD');
    % title('Nyquist rate PSD');
    title('PSD Reconstructed by StOMP');
    grid on;
    subplot(3,1,3);
    plot(f,sig0,'b',f,sig,'r');
    axis([-50 50 -42 0]);
    xlabel('f[MHz]');
    ylabel('PSD');
    title('Comparison');
    grid on;
end
end
figure(4);
N=48:100:M;
plot(N/M,mse(N),'bo-');
axis([0 1 0 0.7]);
title('MSE=E[|| S_{x0} - S_{x} ||_{2}/||S_{x}||_{2}]');
xlabel('Compression rate [N/M]');
ylabel('MSE');
grid on;