bin.rar_direction finding_mayxu6_music_侧向算法_空间谱分析

%MWC检测线性调频信号 clc;clear all;close all; %% 时域测试信号生成 N=5000; % 信号长度200ns fc=900e6; % 信号载频 fs=2.2e9; % 采样频率0.1ns ts=1/fs; % 采样间隔 Ts=1:N; % 采样序列 BW=100e6; BW_K=BW/(N*ts); snr=15; % X0=1*exp(i*(2*pi*fc*Ts*ts)); X1=1*exp(i*(2*pi*fc*Ts*ts+pi*BW_K*ts*ts*Ts.*Ts)); % LFM信号 X=awgn(X1,snr,'measured'); % figure(1); % t=linspace(1,N,N); % plot(t,real(x)); % title('原始信号时域图'); figure(2); f=Ts*fs/N/1e6; plot(f,abs(fft(X))); xlabel('频率(MHz)');ylabel('幅度'); % axis([978 1022 0 4500]); title('原始信号频域图'); %% MWC过程 Mp=20; %伪随机数周期 channel_num=6;% 分支路数 fp=fs/Mp;%伪随机数频率 fs_low=fs/Mp;%压缩采样频率 periodnum=N/Mp;%101 bernoulli_00=zeros(1,2*Mp); xx=zeros(channel_num,N); xxx=zeros(channel_num,N); xx_filter=zeros(channel_num,N); x_sample=zeros(channel_num,N/Mp); % ideal filter ct11=ones(1,periodnum/2); ct111=ones(1,periodnum/2); ct12=zeros(1,N-periodnum); idealfilter=[ct11 ct12 ct111]; bernoulli_1= rand(1,Mp) > 0.5 ; % n是要产生的随机数个数，p是伯努利分布的数学期望 bernoulli_11=ones(1,Mp); bernoulli_11(find(bernoulli_1<0.5))=-1; bernoulli_00=repmat(bernoulli_11,1,2); for k=1:1:channel_num % bernoulli_0= rand(1,Mp) > 0.5 ; % n是要产生的随机数个数，p是伯努利分布的数学期望 % bernoulli0=ones(1,Mp); % bernoulli0(find(bernoulli_0<0.5))=-1; % bernoulli0ct(k,:)=bernoulli0; % bernoulli(k,:)=repmat(bernoulli0,1,periodnum); bernoulli(k,:)=repmat(bernoulli_00(k:k+Mp-1),1,periodnum); xx(k,:)=X.*(bernoulli(k,:));%????????????????? ddd=(fft(xx(k,:))).*idealfilter; xxx(k,:)=ifft(ddd); x_sample(k,:)=downsample( xxx(k,:),Mp); end %% 画MWC的输出 NN=N/Mp; Tss=1:NN; k=3;%某个路 ff=Tss*fs_low/NN; % figure(3); % plot(ff,20*log10(abs(fft(x_sample(k,:)))));xlabel('频率(MHz)');ylabel('幅值(dB)'); % title('输出信号频域图'); % figure(4); % plot((real(x_sample(k,:)))); % title('输出信号时域图'); figure(5); plot(ff,abs(fft(x_sample(k,:))));xlabel('频率(MHz)');ylabel('幅值(dB)'); %% 相位差法求瞬时频率 sig_fai1 = angle(x_sample(k,:)); sig_fai_diff1 = diff(sig_fai1); sig_fai_diff1_adjust = sig_fai_diff1 + (sig_fai_diff1<-pi)*2*pi - (sig_fai_diff1>=pi)*2*pi; win1 = 8; sig_fai_diff1_mean = filter(ones(1,win1)/win1,1,sig_fai_diff1_adjust); sig_inst_fre = sig_fai_diff1_mean/2/pi*fs_low; k1=((sig_inst_fre(99)-sig_inst_fre(90))/10); k2=((sig_inst_fre(89)-sig_inst_fre(80))/10); %% DPT算法 %求调频系数 N1=NN/2; k=1; for n=1:NN-N1 y1(k,n)=x_sample(k,n+N1).*conj(x_sample(k,n)); end x_fft(k,:)=abs(fft(y1(k,:))); [x_max1,tx_max1(k)]=max(x_fft(k,:)); if x_fft(tx_max1-1)<x_fft(tx_max1+1) r1=tx_max1+1; else r1=tx_max1-1; end a1=(r1-tx_max1)*x_fft(r1)/(x_fft(tx_max1)+x_fft(r1)); tx_max=(tx_max1-1+a1)*fs_low/(NN-N1); k_est=tx_max/(N1)*fs_low %求调频系数 figure(6); plot(abs(fft(y1(k,:))));xlabel('频率(MHz)');ylabel('幅值'); %求基带信号初频 for n=1:NN z1(k,n)=exp(-i*(pi*k_est*((1/fs_low).^2)*(n.^2))); end for n=1:NN x_z(k,n)=x_sample(k,n).*z1(k,n); end x_z_fft(k,:)=abs(fft(x_z(k,:))); [x_z_max1,tx_z_max1(k)]=max(x_z_fft(k,:)); if x_z_fft(tx_z_max1-1)<x_z_fft(tx_z_max1+1) r2=tx_z_max1+1; else r2=tx_z_max1-1; end a2=(r2-tx_max1)*x_z_fft(r2)/(x_z_fft(tx_z_max1)+x_z_fft(r2)); tx_z_max=(tx_z_max1-1+a2)*fs_low/(NN); fc_est=tx_z_max-k_est*N1/fs_low/2; figure(7); plot(x_z_fft(k,:));xlabel('频率(MHz)');ylabel('幅值'); %求载频搬移量 for k=1:channel_num X_sample(k,:)=fft(x_sample(k,:)); end ll=zeros(1,channel_num-1); for k=1:channel_num-1 ll(k)=floor(abs(log(X_sample(k+1,:)/X_sample(k,:)))*Mp/(2*pi)); end ll_1=mode(ll(k)) fc_yuan=((ll_1*fs_low+fc_est)) figure(8); plot(fftshift(sig_inst_fre));