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;