基于svd算法实现对海杂波的抑制附Matlab代码

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))),'r-'); %一种转化成分贝的简单方法 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