UWB 超宽带技术 源代码.doc

BIT: 产生原信号二进制比特流PPM-TH: 进行PPM-TH调制Repcode: 进行重复编码TH: 产生TH码TRANSMITTER_2PPM_TH：产生UWB信号waveform: 产生功率归一化的脉冲波形1、FUNCTION: bit% 产生二进制原信号% 原信号比特数numbis作为输入function bits=bit(numbits)bits=rand(1,numbits)0.5;%rand产生的是在01上均匀分布的随机数%这些数0.5的几率各是一半，即bis为0，1的几率各半3、FUNCTION: repcode% 产生重复编码 % Ns ：码元重复数function repbits=repcode(bits,Ns)numbits = length(bits);temprect=ones(1,Ns);temp1=zeros(1,numbits*Ns);temp1(1:Ns:1+Ns*(numbits-1)=bits;temp2=conv(temp1,temprect);repbits=temp2(1:Ns*numbits);3、 FUNCTION: PPM_TH% 引入TH码并进行PPM调制% 参数如下：% seq：二进制源码% fc ：抽样频率% Tc ：时隙长度% Ts ：脉冲平均重复周期% dPPM：PPM引入的时移% THcode ：TH码% 产生两个输出：% 2PPMTHseq ：TH和PPM共同调制信号% THseq ：未经PPM调制的信号function PPMTHseq,THseq = PPM_TH(seq,fc,Tc,Ts,dPPM,THcode)% 调制dt = 1 ./ fc; framesamples = floor(Ts./dt); %每个脉冲的样本数 chipsamples = floor (Tc./dt); PPMsamples = floor (dPPM./dt); THp = length(THcode); totlength = framesamples*length(seq);PPMTHseq=zeros(1,totlength);THseq=zeros(1,totlength);% 引入TH码和PPM %s(t)=sum(p(t-jTs-CjTc-aE)for k = 1 : length(seq) % 脉冲位置,表示第几个脉冲-jTs index = 1 + (k-1)*framesamples; % 引入TH码,-CjTc，表示第几个时隙 kTH = THcode(1+mod(k-1,THp); index = index + kTH*chipsamples; THseq(index) = 1; % 引入PPM时移,-aE，表示在时隙内的位置 index = index + PPMsamples*seq(k); PPMTHseq(index) = 1;end4、 FUNCTION : TH% 产生TH码% Np:跳时码周期% Nh:跳时码最大上界function THcode=TH(Nh,Np);THcode = floor(rand(1,Np).*Nh);5、 FUNCTION : transmitter_2PPM_TH%产生UWB信号% 参数定义如下：% Pow:传输功率% fc;抽样频率% numbits：信号比特数% Ns：每比特脉冲数% Np:跳时码周期% Nh:跳时码最大上界% Ts：脉冲重复周期.% Tc：时隙大小% Tm:脉冲持续时间% tau：脉冲成形因子% dPPM：PPM引入时移% 返回值：% bits:产生比特流% THcode:TH码% Stx：产生信号% ref：未经调制的参照信号function bits,THcode,Stx,ref=transmitter_2PPM_TH(fc,numbit,Ns,Ts,dPPM)% 输入参数Pow = -30; numbits = numbit; Tc = 1e-9; Nh = 10; Np = 5; Tm = 0.5e-9; tau = 0.2e-9; G = 1;% 模拟发射步骤% 二进制原信号bits = bit(numbits);repbits = repcode(bits,Ns); % 重复编码THcode = TH(Nh,Np); % 产生TH码PPMTHseq,THseq = PPM_TH(repbits,fc,Tc,Ts,dPPM,THcode); % 调制% 成形滤波power = (10(Pow/10)/1000; Ex = power * Ts; w0 = waveform(fc,Tm,tau); wtx = w0 .* sqrt(Ex); Sa = conv(PPMTHseq,wtx); Sb = conv(THseq,wtx); % 产生输出信号L = floor(Ts*fc)*Ns*numbits;Stx = Sa(1:L);ref = Sb(1:L);if G %绘图F = figure(1);set(F,Position,30 120 700 400);clatmax = numbits*Ns*Ts;time = linspace(0,tmax,length(Stx);P = plot(time,Stx);set(P,LineWidth,2);ylow=-1.5*abs(min(wtx);yhigh=1.5*max(wtx);axis(0 tmax ylow yhigh);AX=gca;set(AX,FontSize,12);X=xlabel(时间 s);set(X,FontSize,14);Y=ylabel(幅度 V);set(Y,FontSize,14);for j = 1 : numbits tj = (j-1)*Ns*Ts; L1=line(tj tj,ylow yhigh); set(L1,Color,0 0 0,LineStyle, .%间隔比特的线 -,LineWidth,2); for k = 0 : Ns-1 if k 0 tn = tj + k*Nh*Tc; L2=line(tn tn,ylow yhigh); set(L2,Color,0.5 0.5 0.5,LineStyle, .%间隔帧的线 -.,LineWidth,2); end for q = 1 : Nh-1 th = tj + k*Nh*Tc + q*Tc; L3=line(th th,0.8*ylow 0.8*yhigh); set(L3,Color,0.5 0.5 0.5,LineStyle, . :,LineWidth,1); %间隔时隙的线 end endendend6、FUNCTION: waveform% 产生功率归一化的脉冲波形，这里采用的是高斯波形的二阶导数% fc ：抽样频率% Tm ：脉冲持续时间% tau ：成形参数function w0=waveform(fc,Tm,tau);% 产生波形dt = 1 / fc; OVER = floor(Tm/dt); e = mod(OVER,2);kbk = floor(OVER/2);tmp = linspace(dt,Tm/2,kbk);s = (1-4.*pi.*(tmp./tau).2).* . exp(-2.*pi.*(tmp./tau).2);if e % 奇数 for k=1:length(s) y(kbk+1)=1; y(kbk+1+k)=s(k); y(kbk+1-k)=s(k); endelse % 偶数 for k=1:length(s) y(kbk+k)=s(k); y(kbk+1-k)=s(k); endendE = sum(y.2).*dt); w0 = y ./ (E0.5); %功率归一化7、 FUNCTION : Gnoise% 为输入信号引入加性白噪声% 矢量 ebno包含了各个信噪比值% numbits:表示发送端发送的比特数目% 输出:% 经噪声污染的信号output,噪声矩阵noise.function output,noise = Gnoise(input,exno,numpulses)Ex = (1/numpulses)*sum(input.2); %一个单脉冲的平均接收能量ExNo = 10.(exno./10); No = Ex ./ ExNo; nstdv = sqrt(No./2); %噪声的标准差for j = 1 : length(ExNo) noise(j,:) = nstdv(j) .* randn(1,length(input); output(j,:) = noise(j,:) + input;end8、 FUNCTION : IEEEuwb% 根据IEEE 802.15.SG3a.产生信道冲激响应 % fc :抽样频率% TMG:信道增益% 返回:% 1) h0:信道冲激响应% 2) hf:离散信道冲激响应% 3) OT:观察时间OT% 4) ts:离散分辨率% 5) X:信道增益function h0,hf,OT,ts,X = IEEEuwb(fc,ag);% -% 输入参数% -TMG=ag2 % 信道总多径增益OT = 200e-9; % 观测时间 sts = 1e-9; % 离散分辨率 s LAMBDA = 0.0223*1e9; % 簇平均到达因子 (1/s)lambda = 2.5e9; % 簇内脉冲平均到达因子 (1/s)GAMMA = 7.1e-9; % 簇衰减因子 gamma = 4.3e-9; % 簇内脉冲衰减因子sigma1 = 10(3.3941/10); % 簇的信道衰减系数偏差sigma2 = 10(3.3941/10); % 簇内脉冲信道衰减系数偏差sigmax = 10(3/10); % 信道幅度增益的标准偏差% 脉冲衰减阈值，当exp(-t/gamma)rdt时，该脉冲忽略rdt = 0.001;% 峰值阈值 dB，只考虑幅度在峰值-PT范围以内的脉冲PT = 50;G = 1;% -% 簇的形成% -dt = 1 / fc; % 采样频率T = 1 / LAMBDA; % 簇平均到达时间t = 1 / lambda; % 簇内脉冲平均到达时间si = 1;CAT(i)=0; % 第一簇到达时间，初始化为0next = 0; while next OT i = i + 1; next = next + expinv(rand,T); %产生簇的到达时间，服从p(Tn/Tn-1)=lambda*-exp(Tn/Tn-1) if next 0.5)*2-1).*a; ccoeff = sigma1*randn; % 簇衰减 while exp(-next/gamma)rdt pc = pc + 1; next = next + expinv(rand,t); HT(pc) = CT + next; mx = 10*log(omega)-(10*CT/GAMMA)-(10*next/GAMMA); mu = (mx/log(10)-logvar; a = 10(mu+ccoeff+(sigma2*randn)/20); HA(pc) = (rand0.5)*2-1).*a; endend % for i = 1 : NC% 峰值滤波器peak = abs(max(HA);limit = peak/10(PT/10);HA = HA .* (abs(HA)(limit.*ones(1,length(HA);%凡小于limit的脉冲不输出for i = 1 : pc itk = floor(HT(i)/dt); h(itk+1) = HA(i);end% -% 离散相应形式% -N = floor(ts/dt);L = N*ceil(length(h)/N);h0 = zeros(1,L);hf = h0;h0(1:length(h) = h;for i = 1 : (length(h0)/N) tmp = 0; for j = 1 : N tmp = tmp + h0(j+(i-1)*N); end hf(1+(i-1)*N) = tmp;end% 功率归一化E_tot=sum(h.2);h0 = h0 / sqrt(E_tot);E_tot=sum(hf.2);hf = hf / sqrt(E_tot);mux = (10*log(TMG)/log(10) - (sigmax2)*log(10)/20);X = 10(mux+(sigmax*randn)/20);h0 = X.*h0;hf = X.*hf;% -% 图形输出% -if G Tmax = dt*length(h0); time = (0:dt:Tmax-dt); figure(1) S1=stem(time,h0); AX=gca; set(AX,FontSize,14); T=title(信道冲激响应); set(T,FontSize,14); x=xlabel(时间 s); set(x,FontSize,14); y=ylabel(幅度增益); set(y,FontSize,14); figure(2) S2=stairs(time,hf); AX=gca; set(AX,FontSize,14); T=title(离散时间冲激响应); set(T,FontSize,14); x=xlabel(时间 s); set(x,FontSize,14); y=ylabel(幅度增益); set(y,FontSize,14);end9、 FUNCTION: pathloss% 根据给定的距离d，衰减因子gamma以及1米处的信号c0将输入信号衰减% 函数返回衰减后的信号rx以及信道增益attnfunction attn = pathloss(c0,d,gamma)attn = (c0/sqrt(dgamma);10、 FUNCTION : PPMcorrmask_R% 为二进制PPM UWB信号的Rake接收机计算相关掩膜mask% ref：未经PPM调制的参考信号% fc：抽样频率 % numpulses ：传输脉冲数目% dPPM：PPM时移量% rake：离散冲激相应function mask =PMcorrmask_R(ref,fc,numpulses,dPPM,rake)dt = 1 / fc;LR = length(ref);% 功率归一化Epulse = (sum(ref.2).*dt)/numpulses;nref = ref./sqrt(Epulse);% Rake 卷积mref = conv(nref,rake);mref = mref(1:LR);%构造相关掩膜PPMsamples = floor (dPPM ./ dt);sref(1:PPMsamples)=mref(LR-PPMsamples+1:LR);sref(PPMsamples+1:LR)=mref(1:LR-PPMsamples);mask = mref-sref;11、FUNCTION :PPMreceiver% 构造 2PPM TH UWB 信号的接收机，并计算平均错误率BER% R： 表示所使用的波形矩阵，一个波形对应于矩阵的一行% mask：表示相关掩膜 % fc：抽样频率% bits 发射机产生的原始二进制比特流% Ns 每比特的脉冲数（即用几个脉冲表示1比特）% Ts 平均脉冲重复周期，即一帧的长度% 函数返回：% RXbits ：存储经解调后的二进制数据流% BER：存储计算得到的Prb直（误比特率）%function RXbits,BER = PPMreceiver(R,mask,fc,bits,numbit,Ns,Ts)HDSD = 2;% HDSD = 1 - 硬判决，接收机对表示一个比特的Ns个脉冲逐一独立判断。% 将超过门限的脉冲数与低于门限的脉冲数比较，输出较多者对应的比特% HDSD = 2 - 软判决，接收机把Ns个脉冲形成的信号当作一个单独的多脉冲信号% 接收机与信号进行相关的掩膜m(t)是整个比特的脉冲串% 软判决的效果要好于硬判决% 实现相关器% 一个信号波形一行，共有N个波形。每个波形长度为L，即每个波形的样本数numbit=length(bits)N,L = size(R);RXbits = zeros(N,numbit);dt = 1 / fc; framesamples = floor(Ts ./ dt); bitsamples = framesamples * Ns; for n = 1 : N %分别取出N个波形，逐一判断 rx = R(n,:); mx = rx .* mask; if HDSD = 1 % 硬判决 for nb = 1 : numbit mxk = mx(1+(nb-1)*bitsamples:bitsamples+. (nb-1)*bitsamples); No0 = 0; No1 = 0; for np = 1 : Ns mxkp = mxk(1+(np-1)*framesamples:. framesamples+(np-1)*framesamples); zp = sum(mxkp.*dt); %积分周期为Ts if zp 0 No0 = No0 + 1; else No1 = No1 + 1; end end if No0 No1 RXbits(n,nb) = 0; else RXbits(n,nb) = 1; end end % for nb = 1 : numbit end % end of Hard Decision Detection if HDSD = 2 % 软判决 for nb = 1 : numbit mxk = mx(1+(nb-1)*bitsamples:bitsamples+. (nb-1)*bitsamples); zb = sum(mxk.*dt); %积分周期为NsTs if zb 0 % Z0，判断为0 RXbits(n,nb) = 0; else % Z0)&(i0) NF = NF + 1; j = j + 1; index = s_ind(i); % 排序结束后，从尾部，即最大处往前操作 I(j) = index; T(j) = (index-1)*dt; G(j) = hf(index); i = i - 1;end% 计算权重binsamples = floor(ts/dt);if S NF S = NF;endSrake = zeros(1,NF*binsamples);for nf = 1 : NF x = I(nf); y = G(nf); Arake(x) = y; if nf = S Srake(x) = y; % 只输出S个最大的多径分量 end end 13、Function : corrsyn% 将接收信号和参考信号相与，得出时延估计% fc ：抽样频率% 函数返回相与后的结果cfunction C = corrsyn(signal,template,fc)dt = 1 / fc;for s = 1 : length(signal) stmp = signalshift(template,fc,s*dt); % 在所有可能的位置上逐一进行相关运算，将会在原序列的时间偏移处得到最大值。 %因为此时每个谱线都是相互相乘，而没有因为乘以0而受到损失。 C(s) = sum(signal.*stmp);end % for 1 : length(signal)14、Function : signalshift% 在接受参数in中加入移位tfunction out = signalshift(in,fc,t)dt = 1 / fc; % 按给定大小进行时间偏移ss = mod(floor(t/dt),length(in);out = in(end-ss+1:end 1:end-ss);15、 FUNCTION : count_Estimated_Time% 该函数根据估计的时延矩阵计算目标结点到各个参考点的估计距离% 输入参数:% 各个估计时延% 函数返回计算所得1*k的距离矩阵function estimated_Time = count_Estimated_Time(estimated_delay_1,estimated_delay_2,estimated_delay_3);% 构造矩阵err_rangesestimated_Time(1)=estimated_delay_1;estimated_Time(2)=estimated_delay_2;estimated_Time(3)=estimated_delay_3;16、FUNCTION : count_Time% 该函数根据结点距离矩阵计算目标结点到各个参考点的精确时延% 输入参数:% ranges：包含任两点间距离大小为NxN的二维数组% Nx:目标结点（待定位结点）序号% Ref: 长度为k的参考结点数组 % 函数返回计算所得1*k的时延矩阵function Time = count_Time(ranges, Nx, Ref);% 构造矩阵Timek = length(Ref);for i=1:k Time(i) = ranges(Ref(i),Nx)/(3e8);end;17、 FUNCTION : create_network% 在一个正方形区域里产生若干个结点并计算任两点间的距离% 输入参数：% N：结点个数% 正方形边长：area_side% 是否输出图形：G% 输出：% 一个包含各个结点坐标参数的大小为Nx2的二维数组% 一个包含任两点间距离大小为NxN的二维数组function positions, ranges = create_network(N,area_side,G)for i = 1:N positions(i,1)=rand*area_side; positions(i,2)=rand*area_side; j=1; for j=1:(i-1) ranges(i,j)= sqrt(positions(i,1) -. positions(j,1)2 + (positions(i,2) -. positions(j,2)2); while(ranges(i,j)=0) X(i)=rand*50; Y(i)=rand*50; ranges(i,j)=sqrt(positions(i,1)-. positions(j,1)2+(positions(i,2)-. positions(j,2)2); end ranges(j,i)=ranges(i,j); endendif G scatter(positions(:,1),positions(:,2),20,filled); axis(0 area_side 0 area_side); xlabel(X m); ylabel(Y m); box on;end18、FUNCTION : find_position% 根据测距估计值定位% 输入参数:% positions：包含各个结点坐标参数的大小为Nx2的二维数组% ranges：包含任两点间距离大小为NxN的二维数组% Nx:目标结点（待定位结点）序号% Ref: 长度为k的参考结点数组 % sigma_2:测距误差% G:是否输出图形% 函数返回计算所得的目标结点位置并参照其实际位置计算定位误差% function PosNx, ErrNx = find_LSE_position(positions, ranges, Nx, Ref,estimated_Time,G);% 用估计的时延k = length(Ref);for i=1:kerr_ranges(i)=estimated_Time(i)*3e8;end% 构造矩阵Afor i=1:(k-1) A(i,1) = positions(Ref(i),1) - positions(Ref(k),1); A(i,2) = positions(Ref(i),2) - positions(Ref(k),2);endA=-2*A;% 构造矩阵bb=zeros(2,1);for i=1:(k-1) b(i) = err_ranges(i)2 -. err_ranges(k)2 - positions(Ref(i),1)2 +. positions(Ref(k),1)2 - positions(Ref(i),2)2 +. positions(Ref(k),2)2;end% 计算目标结点位置PosNx=Ab;% 计算定位误差ErrNx = sqrt(PosNx(1)-positions(Nx,1)2+(PosNx(2)-. positions(Nx,2)2);% 图形输出if G scatter(positions(:,1),positions(:,2); xlabel(X m); ylabel(Y m); box on; hold on; scatter(PosNx(1), PosNx(2), 200, filled, k,p); scatter(positions(Nx,1),positions(Nx,2),200,. filled,); for i=1:k scatter(positions(Ref(i),1),positions(Ref(i),2),. filled,r,s); end hold off;end19、 FUNCTION : find_LSE_position% 根据节点间距精确定位% 输入参数:% positions：包含各个结点坐标参数的大小为Nx2的二维数组% ranges：包含任两点间距离大小为NxN的二维数组% Nx:目标结点（待定位结点）序号% Ref: 长度为k的参考结点数组 % sigma_2:测距误差% G:是否输出图形% 函数返回计算所得的目标结点位置并参照其实际位置计算定位误差function PosNx, ErrNx = find_LSE_position(positions, ranges, Nx, Ref,sigma_2, G);% 加入误差N = size(ranges,1);err_ranges = ranges + sqrt(sigma_2)*randn(N);% 构造矩阵Ak = length(Ref);for i=1:(k-1) A(i,1) = positions(Ref(i),1) - positions(Ref(k),1); A(i,2) = positions(Ref(i),2) - positions(Ref(k),2);endA=-2*A;% 构造矩阵bb=zeros(2,1);for i=1:(k-1) b(i) = err_ranges(Ref(i),Nx)2 -. err_ranges(Ref(k),Nx)2 - positions(Ref(i),1)2 +. positions(Ref(k),1)2 - positions(Ref(i),2)2 +. positions(Ref(k),2)2;end% 计算目标结点位置