function y_new = SVD(num,L)
% 对某一距离单元的回波进行预处理，并返回预处理后的数据
%%% num表示要取第几个距离单元的时域数据 %%%
%%% L表示Hankel矩阵的列 %%%
%%% y_new表示与处理后的时域数据 %%%
%% 导入并测试需要处理的数据
%小数据分析
load Time_dataS.mat;
%T = 0.149;
f = 3700000; % 高频地波雷达的频率
T = 1/f;
Num_T =1024;
doppler = 1/T/Num_T*(-Num_T/2:Num_T/2-1);
t = 1:1:1024;
y = Yt(30,:);%data30
yf = db(abs(fftshift(fft(y))));
%大数据分析
% load data75V8;
% t = 1:1:16384;
% FileNumber = 8;
% sf = db(abs(fftshift(fft(s))));
% f = 3700000;
% T = 1/f;
% Num_T = 2048*FileNumber;
% fd = 1/T/Num_T*(-Num_T/2:Num_T/2-1);
% figure(1);
% plot(t,s);title('时域数据');
% figure(2);
% plot(fd,sf);title('频域数据');axis([-0.2e5,0.2e5,-170,-90]);
%% 奇异值分解处理
%% 构造汉克尔矩阵"H"
L = 60; % L是输出值，也就是汉克尔矩阵的列（L取125效果最佳）
N = length(y); % 时域数据长度
R = N-L+1; % 汉克尔矩阵的行
H = [];
for i = 1:R
for j = 1:L
H(i,j) = y(i+j-1);
end
end
%% 对汉克尔矩阵进行奇异值分解，并对相应奇异值置零
[U,S,V] = svd(H); % 对Hankel矩阵进行奇异值分解，S中包含奇异值矩阵
s1 = S;
s1(1,1) = 0; % 对应的奇异值置零
for i = 2:5
s1(i,i) = 0;
end
sigma=diag(S); %获得主对角线元素
figure(1);stem(sigma);
title('奇异值分布');
for i = 8:L
s1(i,i) = 0;
end
u1 = U; u1(1,1) = 0; % U,V也要相应的变化
v1 = V; v1(1,1) = 0;
for i = 2:5
u1(i,i) = 0;
v1(i,i) = 0;
end
for i = 8:L
u1(i,i) = 0;
v1(i,i) = 0;
end
H_new = u1*s1*v1'; % 构造新的汉克尔矩阵
%% 根据新的汉克尔矩阵重新构造时域数据
p = 0;
for k = 1:N
y_new(k) = 0;
if k<L
cl = k;
for ro = 1:1
y_new(k) = H_new(ro,cl) + y_new(k);
cl = cl-1; %控制列数
end
y_new(k) = y_new(k)/k; %将第一行前L列H_new元素分别/(1:L)例H_new(1,2)/2 赋值给 y_new共L-1个数
elseif k >= L & k <= N-L+1
for ro = 1:L
y_new(k) = H_new(ro+p,L-ro+1)+y_new(k);%将L;N-L+1个数赋值，若k=L,将H_new前L行L列数相加
end
p = p+1;
y_new(k) = y_new(k)/L;
elseif k > N-L+1
cl = L;
for ro = k-L+1:N-L+1
y_new(k) = H_new(ro,cl)+y_new(k);
cl = cl-1;
end
y_new(k) = y_new(k)/(N-L+1);
end
end
%% 作图
N1 = 1024; %number of zero padding
figure(2);
t = 1:N1;
plot(t,y);title('预处理前的时域数据');xlabel('t');ylabel('海杂波数据')
Y = fftshift(fft(y,N1));%Y=fft(y,N1);Y=fftshift(Y);
Y_new = fftshift(fft(y_new,N1));
figure(3);
plot(t,y_new);title('预处理后的时域数据');xlabel('t');ylabel('海杂波数据');
figure(4);
plot(doppler,db(abs(Y)/max(abs(Y))),'--'); %一种转化成分贝的简单方法
hold on;
plot(doppler,db(abs(Y_new)/max(abs(Y_new))));axis([-0.5e6,0.5e6,-80,0]);
% figure(5);
% plot(doppler,db(abs(Y)/max(abs(Y))),'--'); %一种转化成分贝的简单方法
% xlabel('频率');ylabel('幅度');axis([-0.5e6,0.5e6,-80,0]);
title('干扰信号抑制前后频谱对比');
legend('干扰信号抑制前','干扰信号抑制后')
xlabel('f (Hz)'),ylabel('Am (dB)')
end
评论7