www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/simulation_sqy4Hilbert.m
clc
clear
close all
%%%******** 读数据 ****************%%%%
%%------------------------------------------------------------------------
input_file_name_real = 'e:\realD74_test.dat';
input_file_name_img = 'e:\imagD74_test.dat';
output_file_name_real = 'e:\0917data_real_D74.dat';
output_file_name_img = 'e:\0917data_img_D74.dat';
% output_file_name_bp = 'E:\GPS\data_bp.dat';
%%------------------------read----------------------------------------------
read_fid_real = fopen( input_file_name_real, 'r');
if read_fid_real ==-1
disp( 'can not open the file1 for read !, please check the file name and path !' );
end
read_fid_img = fopen( input_file_name_img, 'r');
if read_fid_img ==-1
disp( 'can not open the file2 for read !, please check the file name and path !' );
end
%%------------------------write-------------------------------------------
% write the data to file
write_fid_real = fopen(output_file_name_real, 'w' );
if write_fid_real == -1
disp('can not open the file for writing !, please check the file name and path !' );
end
write_fid_img = fopen(output_file_name_img, 'w' );
if write_fid_img ==-1
disp('can not open the file2 for writing !, please check the file name and path !' );
end
% write_fid_bp = fopen(output_file_name_bp, 'w' );
% if write_fid_bp ==-1
% disp('can not open the file2 for writing !, please check the file name and path !' );
% end
%%-------------------------------------------------------------------------
N = 400;
count = 60000;
input_buffer_real = zeros(1, count);
input_buffer_img = zeros(1, count);
%%-------------------自适应-------------------------------------------------
ii=0;
BlockNo=80;
while (BlockNo~=0)
for i=1:4 %读数据60000*2*4
[input_buffer_real, count] = fread(read_fid_real, count, 'single');
[input_buffer_img, count] = fread(read_fid_img, count, 'single');
Real_channel(i,:) = input_buffer_real.'; %+j*input_buffer2';
Img_channel(i,:) = input_buffer_img.';
% Img_channel(i,:)=imag(hilbert(input_buffer_real));
% Data_channel(i,:) = input_buffer_real.'+j*input_buffer_img.';
end
% [write_count_bp] = fwrite(write_fid_bp, Real_channel(1,:), 'int8');
if BlockNo==80
% for i=1:count/N
for i =1:1
for kx=1:4
xxr(kx,:) = Real_channel(kx,1+(i-1)*N:i*N);% 每次处理400个点
xxi(kx,:) = Img_channel(kx,1+(i-1)*N:i*N);
end
%%%%%************ 求相关阵 ***************
Rz_real=zeros(4,4);
Rz_img=zeros(4,4);
for k=1:N
for m=1:4
for n=m:4
Rz_real(m,n) = Rz_real(m,n)+(xxr(m,k)*xxr(n,k)+xxi(m,k)*xxi(n,k));
Rz_img(m,n) = Rz_img(m,n)-(xxr(m,k)*xxi(n,k)-xxr(n,k)*xxi(m,k));
end
end
end
for m=2:4
for n=1:m-1
Rz_real(m,n) = Rz_real(n,m); %%%%实对称
Rz_img(m,n) = -Rz_img(n,m);
end
end
%%%%%%%******************* 嵌入DSP **********************************%%%%%%
%%%%%%%%************* 对角加载 **************************
% % diagload=sqrt(trace(Rz_real)/rank(Rz_real))*0.01;
% % Rz_real=Rz_real+diagload*eye(size(Rz_real,2));
tra = Rz_real(1,1)+Rz_real(2,2)+Rz_real(3,3)+Rz_real(4,4);
for mm = 1:4
Rz_real(mm,mm) = Rz_real(mm,mm) + tra/100;
end
% Rz = Rz_real + j*Rz_img ; %%%求相关矩阵,精度差小于1,扩大了2的16倍
%%%%************** Cholesky分解************%%%
cr=zeros(4);
lr=zeros(4);
ci=zeros(4);
li=zeros(4);
for k=1:4
cr(k,1)=Rz_real(k,1);
ci(k,1)=Rz_img(k,1);
end
dr1=cr(1,1);
for k=1:4
lr(k,1)=cr(k,1)/dr1; %扩大Nshift倍
li(k,1)=ci(k,1)/dr1;
end
for k=2:4
cr(k,2)=Rz_real(k,2)-(cr(k,1)*lr(2,1)+ci(k,1)*li(2,1));
ci(k,2)=Rz_img(k,2)+(cr(k,1)*li(2,1)-ci(k,1)*lr(2,1));
end
dr2=cr(2,2);
for k=2:4
lr(k,2)=cr(k,2)/dr2; %扩大Nshift倍
li(k,2)=ci(k,2)/dr2;
end
for k=3:4
cr(k,3)=Rz_real(k,3)-(cr(k,1)*lr(3,1)+ci(k,1)*li(3,1)+cr(k,2)*lr(3,2)+ci(k,2)*li(3,2));
ci(k,3)=Rz_img(k,3)+(cr(k,1)*li(3,1)-ci(k,1)*lr(3,1)+cr(k,2)*li(3,2)-ci(k,2)*lr(3,2));
end
dr3=cr(3,3);
for k=3:4
lr(k,3)=cr(k,3)/dr3;
li(k,3)=ci(k,3)/dr3;
end
cr(4,4)=Rz_real(4,4)-(cr(4,1)*lr(4,1)+ci(4,1)*li(4,1)+cr(4,2)*lr(4,2)+ci(4,2)*li(4,2)+cr(4,3)*lr(4,3)+ci(4,3)*li(4,3));
ci(4,4)=Rz_img(4,4)+(cr(4,1)*li(4,1)-ci(4,1)*lr(4,1)+cr(4,2)*li(4,2)-ci(4,2)*lr(4,2)+cr(4,3)*li(4,3)-ci(4,3)*lr(4,3));
lr(4,4)=1;
li(4,4)=0;
% c1=cr+j*ci; %%%c缩小Nshift
% l1=lr+j*li; %%%l扩大Nshift
%%% ------------------------定点求Y ----------------------------
yr=zeros(4,1);
yi=zeros(4,1);
yr(1,1)= 1/cr(1,1); %%扩大Nshift
yi(1,1)=0;
yr(2,1)= (-cr(2,1)*yr(1,1)+ci(2,1)*yi(1,1)) /cr(2,2);
yi(2,1)= (-cr(2,1)*yi(1,1)-ci(2,1)*yr(1,1))/cr(2,2);
yr(3,1)= (-cr(3,1)*yr(1,1)+ci(3,1)*yi(1,1)-cr(3,2)*yr(2,1)+ci(3,2)*yi(2,1))/cr(3,3);
yi(3,1)= (-cr(3,1)*yi(1,1)-ci(3,1)*yr(1,1)-cr(3,2)*yi(2,1)-ci(3,2)*yr(2,1))/cr(3,3);
yr(4,1)= (-cr(4,1)*yr(1,1)+ci(4,1)*yi(1,1)-cr(4,2)*yr(2,1)+ci(4,2)*yi(2,1)-cr(4,3)*yr(3,1)+ci(4,3)*yi(3,1))/cr(4,4);
yi(4,1)= (-cr(4,1)*yi(1,1)-ci(4,1)*yr(1,1)-cr(4,2)*yi(2,1)-ci(4,2)*yr(2,1)-cr(4,3)*yi(3,1)-ci(4,3)*yr(3,1))/cr(4,4);
% y1=yr+j*yi;
%%% ------------------------定点求W ----------------------------
wr=zeros(4,1);
wi=zeros(4,1);
wr(4,1)=yr(4,1);
wi(4,1)=yi(4,1);
wr(3,1)=yr(3,1)-(lr(4,3)*wr(4,1)+li(4,3)*wi(4,1));
wi(3,1)=yi(3,1)-(lr(4,3)*wi(4,1)-li(4,3)*wr(4,1));
wr(2,1)=yr(2,1)-(lr(4,2)*wr(4,1)+li(4,2)*wi(4,1)+lr(3,2)*wr(3,1)+li(3,2)*wi(3,1));
wi(2,1)=yi(2,1)-(lr(4,2)*wi(4,1)-li(4,2)*wr(4,1)+lr(3,2)*wi(3,1)-li(3,2)*wr(3,1));
wr(1,1)=yr(1,1)-(lr(4,1)*wr(4,1)+li(4,1)*wi(4,1)+lr(3,1)*wr(3,1)+li(3,1)*wi(3,1)+lr(2,1)*wr(2,1)+li(2,1)*wi(2,1));
wi(1,1)=yi(1,1)-(lr(4,1)*wi(4,1)-li(4 ,1)*wr(4,1)+lr(3,1)*wi(3,1)-li(3,1)*wr(3,1)+lr(2,1)*wi(2,1)-li(2,1)*wr(2,1));
% wr=wr/wr(1,1);
% wi=wi/wr(1,1);
%%%----------非连续自适应-----------------------------------------------------------
end
end
output_buffer_real = (real((wr+j*wi)'*(Real_channel+j*Img_channel))>0)*2-1;
output_buffer_img = (imag((wr+j*wi)'*(Real_channel+j*Img_channel))>0)*2-1;
[write_count_real] = fwrite(write_fid_real, output_buffer_real, 'int8');
[write_count_img] = fwrite(write_fid_img, output_buffer_img, 'int8');
%%%%------------连续自适应-------------------------------------------
% output_buffer_real = floor(real((wr+j*wi)'*(xxr+j*xxi)*2^21));
% output_buffer_img = floor(imag((wr+j*wi)'*(xxr+j*xxi)*2^21));
% [write_count_real] = fwrite(write_fid_real, output_buffer_real, 'int8');
% [write_count_img] = fwrite(write_fid_img, output_buffer_img, 'int8');
% end
%%%-----------------------------------------------------------------------
BlockNo=BlockNo-1
end
fclose(write_fid_real);
fclose(write_fid_img);
fid1 = fopen('D:\sample\data_bp.dat');
bpro_data = fread(fid1,60000).';
fid2 = fopen('D:\sample\data_img_400.dat');
pro_data = fread(fid2,60000).';
Q = 60000;
fs = 6e6;
figure(1)
subplot(2,1,1),plot((-Q/2:Q/2-1)/Q*fs,10*log10(fftshift((abs(fft(bpro_data))))))
axis([-3e6,3e6,0,80])
title('before processing')
subplot(2,1,2),plot((-Q/2:Q/2-1)/Q*fs,10*log10(fftshift((abs(fft(pro_data))))))
axis([-3e6,3e6,0,80])
title('after processing(200)')