clc;clear
B = 15.0e6;
T = 10.e-6;
Fs=2*B;
f0=8e7;
mu = B / T;
N=1024;
t = linspace(-T/2., T/2., N);
LFM=exp(j*2*pi*(mu .* t.^2 / 2.));
LFMFFT = fftshift(fft(LFM));
%freqlimit = 0.5 / 1.e-9;
%freq = linspace(-freqlimit/1.e6,freqlimit/1.e6,N);
figure(1)
plot(t*1e6,LFM,'k');
axis([-1 1 -1.5 1.5])
grid;
xlabel('时间/us')
ylabel('幅度/v')
n=1:N;
ts=n/Fs;
St=LFM;
figure(1)
subplot(211)
plot(t*1e6,St);
xlabel('Time in u sec');
title('线性调频信号');
grid on;axis tight;
subplot(212)
freq=linspace(-Fs/2,Fs/2,N);
plot(freq*1e-6,fftshift(abs(fft(St))));
xlabel('Frequency in MHz');
title('线性调频信号的幅频特性');
grid on;axis tight;
z=fft(St,N);
figure(2)
subplot(211)
freq=linspace(-Fs/2,Fs/2,N);
mag=abs(z);
f=(0:(length(z)-1))*Fs/length(z);
plot(f,mag);
xlabel('Frequency in MHz');
title('线性调频信号的幅频特性');
grid on;axis tight;
K=100; % 稀疏度(做FFT可以看出来)
M=512; % 测量数(M>=K*log(N/K),至少40,但有出错的概率)
x=St;
%% 2. 时域信号压缩传感
Phi=randn(M,N); % 测量矩阵(高斯分布白噪声)
s=Phi*x.'; % 获得线性测量
zs=fft(s,N);
figure(6)
freq=linspace(-Fs/2,Fs/2,N);
mags=abs(zs);
ff=(0:(length(zs)-1))*Fs/length(zs);
plot(ff,mags);
xlabel('Frequency in MHz');
title('线性调频信号的幅频特性');
grid on;
axis tight;
%% 3. 正交匹配追踪法重构信号(本质上是1-范数最优化问题)
m=2*K; % 算法迭代次数(m>=K)
Psi=fft(eye(N))/sqrt(N); % 傅里叶正变换矩阵
T=Phi*Psi'; % 恢复矩阵(测量矩阵*正交反变换矩阵)
hat_y=zeros(1,N); % 待重构的谱域(变换域)向量
Aug_t=[]; % 增量矩阵(初始值为空矩阵)
r_n=s; % 残差值
for times=1:m; % 迭代次数(有噪声的情况下,该迭代次数为K)
for col=1:N; % 恢复矩阵的所有列向量
product(col)=abs(T(:,col)'*r_n); % 恢复矩阵的列向量和残差的投影系数(内积值)
end
[val,pos]=max(product); % 最大投影系数对应的位置
Aug_t=[Aug_t,T(:,pos)]; % 矩阵扩充
T(:,pos)=zeros(M,1); % 选中的列置零（实质上应该去掉，为了简单我把它置零）
aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; % 最小二乘,使残差最小
r_n=s-Aug_t*aug_y; % 残差
pos_array(times)=pos; % 纪录最大投影系数的位置
end
hat_y(pos_array)=aug_y; % 重构的谱域向量
figure(2)
subplot(212)
f=(0:(length(hat_y)-1))*Fs/length(hat_y);
plot(f,(abs(hat_y)));
title('CS result of using dechirp-basis');
hat_x=real(Psi'*hat_y.'); % 做逆傅里叶变换重构得到时域信号
%% 4.恢复信号和原始信号对比
figure(3)
%hold on;
subplot(311)
plot(hat_x,'b') % 重建信号
subplot(212)
plot(x,'r')
figure(4)
subplot(311)
plot(t*1e6,hat_x,'b'); % 重建信号
%axis([-5 5 -1 1])
xlabel('Time in u sec');
title('重建信号');
legend('Recovery')
grid on;
subplot(312)
plot(t*1e6,x,'r');
%plot(t*1e6,x); % 原始信号
xlabel('Time in u sec');
title('原始信号');
legend('Original')
grid on;axis tight;
xx=x-hat_x';
subplot(313)
plot(t*1e6,xx);
title('信号误差');
%norm(hat_x.'-x)/norm(x) % 重构误差
