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

    clc;
clear all;
close all;
sig01=load('shuju3.txt');
sig0=sig01;
xfs=50;
points=2048;
f = xfs*(-points:points-1)/points;

figure;
plot(f,sig0);
grid on;
M = 2^11;
N = 900;

%sig0 = MakeBlocks(M)';
qmf = MakeONFilter('Haar',1);

SNR =30;
y=awgn(sig0,SNR);

A = MatrixEnsemble(N,M);
AM=MatrixEnsemble(2*N,2*M);
L =3;
alpha0 = FWT_PO(y, L, qmf);

% Generate Random Dictionary
% Phi = MatrixEnsemble(N,M);
Phi=AM;

% generate the vector S
S = Phi * alpha0;

% Solve the BP problem
alpha = SolveMP(Phi, S,2*M);%%SolveOMP,SolveMP,SolveBP
sig = IWT_PO(alpha,L,qmf);

for ii=1:length(sig0)
%     MSE(ii)=(sig0(ii)-sig(ii))^2/sig0(ii)^2;
      MSE(ii)=abs(sig0(ii)-sig(ii));
end

figure;
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;

% figure;
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 BP');
grid on;
subplot(3,1,3);
plot(f,sig0,'b',f,sig,'r');
axis([-50 50 -42 0]);
xlabel('f[MHz]');
ylabel('PSD');
title('Compare');
grid on;