pswf_numerical.m

用特征值分解法获得长椭球波脉冲

%#####################################################
%% 利用特征值分解算法获得PSWF
%  矩阵H的特征向量v(:,i)表示H的特征值d(i)对应的PSWF
%% 参考文献：A Novel Ultra-Wideband Design Algorithm
%          IEEE Communication Letters, 2003.5
%% Author： CHI YG， WF
%% Date：2006.12.15正确的！！！该目下的figure文件由本程序产生
%#####################################################

clear;clc;close all;
%hold on;

%%###程序需要设定的参数
a=3.1e+9;
b=10.6e+9;           %设定H(f)的频带范围[a,b]
t_half=0.5e-9;       %PSWF的时宽，需要凭经验确定，或根据程序得到的波形来调整
N_samp=129;           %矩阵H：N_samp*N_samp
t_samp=2*t_half/(N_samp-1);         %保证采样间隔为1/t_samp
%##################

%以下完成最初时域函数表达式绘图工作
t=-t_half:1e-12:t_half; %设定t值范围以及采样间隔，为了画的效果，采样较小
y=2*b*sinc(2*b*t)-2*a*sinc(2*a*t);%生成滤波器的时域函数表达式

%以下完成对时域表达式64抽样的工作
N=2*N_samp-1;        %为了得到N_samp阶的矩阵H，需采样N点
x=linspace(-2*t_half,2*t_half,N);  %在以上t值区间对t采样129个点
h=zeros(1,N); %生成一个1行N列矩阵，用来存储64抽样的结果
for n=1:N;  %设定n的循环区间
h(n)=2*b*sinc(2*b*x(n))-2*a*sinc(2*a*x(n));%获得采样点值
end

% 以下完成最后要分析的离散化矩阵H的生成工作
H=zeros(N_samp,N_samp);     %生成一个65行65列的全零矩阵，用来存储离散化以后的矩阵
for i=1:N_samp;
    for j=1:N_samp;
        H(i,j)=h(N_samp+abs(i-j));
        H(j,i)=h(N_samp+abs(i-j));     %以上循环完成64点抽样值到矩阵的转换
    end 
end  %察看生成的离散化矩阵
       
%以下完成最后的特征向量绘图工作
[v,d]=eigs(H,6);        %获得其6个最大的特征值d及其特征向量v
y=linspace(-t_half,t_half,N_samp);
figure;
subplot(2,2,1);plot(y*1e9,v(:,1));grid;xlabel('Time (ns)','FontSize',12);ylabel('Amplitude (v)','FontSize',12);   %绘制出特征向量1的时域波形图
subplot(2,2,2);plot(y*1e9,v(:,2));grid;xlabel('Time (ns)','FontSize',12);ylabel('Amplitude (v)','FontSize',12);  %绘制出特征向量2的时域波形图
subplot(2,2,3);plot(y*1e9,v(:,3));grid;xlabel('Time (ns)','FontSize',12);ylabel('Amplitude (v)','FontSize',12);
subplot(2,2,4);plot(y*1e9,v(:,4));grid;xlabel('Time (ns)','FontSize',12);ylabel('Amplitude (v)','FontSize',12);


%The following are Autocorrelation and Cross-correlation
g1=(v(:,1))'; g2=(v(:,2))';g3=(v(:,3))'; 
g4=(v(:,4))';
rg1=conv(g1,g1);     %g1的自相关
rg2=conv(g1,g2);     %g1、g2的互相关
z=linspace(-1e-9,1e-9,length(rg1));
figure;
subplot(1,2,1);plot(z*1e9,rg1);grid ;xlabel('Time (ns)','FontSize',12);ylabel('Amplitude (v)','FontSize',12);
subplot(1,2,2);plot(z*1e9,rg2);grid ;xlabel('Time (ns)','FontSize',12);ylabel('Amplitude (v)','FontSize',12);


%The following is the power spectral
%######## FFT的点数改变时，需要改变m ################
m=linspace(-8e-9,8e-9,2048); % 保证采样间隔为t_samp
[AW1,f]=cftbyfft(g1,m);Aw1=20*log10(AW1/max(abs(AW1))); 
[AW2,f]=cftbyfft(g2,m);Aw2=20*log10(AW2/max(abs(AW2)));
[AW3,f]=cftbyfft(g3,m);Aw3=20*log10(AW3/max(abs(AW3))); 
[AW4,f]=cftbyfft(g4,m);Aw4=20*log10(AW4/max(abs(AW4)));
[f1,EIRP_indoor]=FCCmask;          %FCC对室内UWB系统发射功率谱的规定
EIRP_indoor=EIRP_indoor+41.3;
figure;
subplot(2,1,1);plot (f*1e-9,Aw1,'--k'); %画出数值解法得到的第一个PSWF脉冲
axis([1.5 12.5 -120 0]);
grid;xlabel('Frequency (GHz)','FontSize',12); ylabel('|P(f)|^2 (dBw/Hz)','FontSize',12);

subplot(2,1,2);plot (f1,EIRP_indoor,'k'); %画出EIRP_indoor
axis([1.5 12.5 -120 -40]);
grid;xlabel('Frequency (GHz)','FontSize',12); ylabel('|P(f)|^2 (dBw/Hz)','FontSize',12);

figure;
subplot(2,2,1); plot(f*1e-9,Aw1);
axis([1.5 12.5 -120 0]);grid;xlabel('Frequency (GHz)','FontSize',12); ylabel('|P(f)|^2 (dBw/Hz)','FontSize',12);
subplot(2,2,2); plot(f*1e-9,Aw2);
axis([1.5 12.5 -120 0]);grid;xlabel('Frequency (GHz)','FontSize',12); ylabel('|P(f)|^2 (dBw/Hz)','FontSize',12);
subplot(2,2,3); plot(f*1e-9,Aw3);
axis([1.5 12.5 -120 0]);grid;xlabel('Frequency (GHz)','FontSize',12); ylabel('|P(f)|^2 (dBw/Hz)','FontSize',12);
subplot(2,2,4); plot(f*1e-9,Aw4);
axis([1.5 12.5 -120 0]);grid;xlabel('Frequency (GHz)','FontSize',12); ylabel('|P(f)|^2 (dBw/Hz)','FontSize',12);

%####### 统计能量集中程度energy_efficiency>99%时脉冲波形的时宽tp########
for i=1:N_samp
    g11(i)=abs(g1(i))^2;   %第一个波形
    %g11(i)=abs(g2(i))^2;  %第二个波形
    %g11(i)=abs(g3(i))^2;
    %g11(i)=abs(g4(i))^2; 
    %g11(i)=abs(g5(i))^2;
    %g11(i)=abs(g6(i))^2;
end
sum1=sum(g11);
energy_efficiency=0;
i=(N_samp+1)/2+1;
sum11=g11((N_samp+1)/2);
while energy_efficiency<0.99
    sum11=sum11+2*g11(i);
    energy_efficiency=sum11/sum1;
    i=i+1;
end
temp=i;
tp=2*(y(temp));%能量集中系数〉99%对应的脉冲时宽。
B_computation=1/tp;