% load 000time_signal.txt
%%物理参数
pi=3.1415926;
time_AUtoSI = 2.418884326505e-17; %时间单位转换
time_SItoAU = 1.0/time_AUtoSI; %时间单位转换
fre_AUtoSI = 1/time_AUtoSI; %频率单位转换
fre_SItoAU = time_AUtoSI;
%% ========================= 激光参数 ========================= %%
total=15.99; % 总的光周期数目
aLamda1=346; % 波长
WindowWidth=1.2; % 窗大小
N_Shown_Timepoints=151; % 时频图显示的格点数
t_Initial=0; % 时频图显示的时间开始范围
t_Final=6; % 时频图显示的时间结束范围
HarmonicInitial=3; % 时频图显示频率开始范围
HarmonicFinal=25; % 时频图显示频率结束范围
HarmonicGiven=15;
%% ========================= 激光参数 ========================= %%
OMGA1=2*pi*2.99792458*100/aLamda1*1e15*fre_SItoAU; %基频光圆频率
load t_at.txt; % 读入数据
a1=t_at(:,1); % 1列时间
a2=t_at(:,2); % 2列时阈电流
%nd = size(a1);
nline=length(a1); % 数组长度
ntw=nline
OMGA_Domain = zeros(1,ntw); % 由时间格点数组 => 频率格点数组
%% ========================= 时域电流 ========================= %%
figure;
plot(a1(1:length(a1))/2,a2(1:length(a2)) ); % 时域电流
%% ========================= interband ========================= %%
dt=a1(2)-a1(1); % 单位时间格点
% dt is originally in units of o.c.; therefore, we need to transform to
% a.u.
dti_AU=dt/OMGA1*2*pi; % 单位时间格点 （o.c.=>原子单位制）
a1=a1/OMGA1*2*pi;
% maximal frequency
fmax = pi/dti_AU; % 单位时间格点=>最大频率
a2_OMGA = fft(a2); % 时域电流=>fft到频域电流
ntw=nline;
%% ========================= 时频分析时间格点 ========================= %% 时频分析只要部分时域信号
Window= 1.0/(WindowWidth*OMGA1);
% It_Initial = floor((t_Initial*2*pi -(-total/2.0)*2*pi)/OMGA1/dti_AU)+1 % 时频分析开始时间格点
% It_Final = floor((t_Final*2*pi -(-total/2.0)*2*pi)/OMGA1/dti_AU) +1 % 时频分析结束时间格点
It_Initial = floor((t_Initial*2*pi )/OMGA1/dti_AU)+1 % 时频分析开始时间格点
It_Final = floor((t_Final*2*pi )/OMGA1/dti_AU) +1 % 时频分析结束时间格点
%% Grids Number in time domain
ntioc = floor((ntw)/total+0.51); % 每个周期内时间格点数目
N_Shown_Timepoints_Range = floor(abs((t_Initial-t_Final)*ntioc)+0.51) % 时频分析用格点数目<小于真实时间格点
N_DeltaTi=floor((N_Shown_Timepoints_Range)/floor(N_Shown_Timepoints)+0.51) % 每隔几个真实时间格点做傅里叶变换
%% ========================= 时频分析频率格点 =========================
N_Shown_Energy_Range = N_Shown_Timepoints;
HarmonicFinal=HarmonicFinal*OMGA1; HarmonicInitial=HarmonicInitial*OMGA1;
DeltaHarmonic = abs(HarmonicFinal-HarmonicInitial)/(N_Shown_Energy_Range)
%% ========================= 时频分析频率格点 =========================
%% ========================= 设置频率格点 =========================
dwi=(2.0*fmax)/(ntw); % 频率间隔
OMGA_Domain(1:ntw/2+1)=(0:ntw/2)*dwi; % 频率格点
OMGA_Domain(ntw/2+2:ntw)=dwi*(1:(ntw)/2-1)-fmax;
%% ========================= interband ========================= %%
figure;
plot(a1(1:length(a1))/2,a2(1:length(a2)) );
%% ========================= 频谱 ========================= %%
figure(2)
Harm_OMGA=abs(a2_OMGA).*abs(a2_OMGA); % 频谱
semilogy(abs(OMGA_Domain(floor(2):floor(0.5*ntw))/OMGA1),abs(Harm_OMGA(floor(2):floor(0.5*ntw))),'g' )
set(gca,'Position',[.12 .17 .80 .74]);
set(gcf,'position',[100 100 580 460],'Units','centimeters');
axis([0 33, 1e-23 1e10]);
set(gca,'xtick',[1 9 17 25 33]);set(gca,'ytick',[1e-20 1e-15 1e-10 1e-5 1 1e5]);
xlabel('\fontsize{16}High harmonic order');
ylabel('\fontsize{16}Spectral intensity(a.u.)')
set(gca,'fontsize',12); set(gca,'FontName','Helvetica');
set(gca,'linewidth',2); set(get(gca,'Children'),'linewidth',2);
%title('\fontsize{20} inter&intra')
saveas(gcf,'harmonic_yield.jpg')
%close(2)
%% ========================= 小窗时频分析 ========================= %%
AtwInteger=0.0;
fid1 = fopen('time_omga_yieldlog.txt','wt'); %log-scale files
fid2 = fopen('time_omga_yield.txt','wt'); %linear-scale files
for ie=1:N_Shown_Energy_Range
HARMONICOMGA = HarmonicInitial+ie*DeltaHarmonic;
if rem(ie,10) == 1; disp(['Number of steps: ' , num2str(ie),' / ', num2str(N_Shown_Energy_Range)]); end
nj2=0;
for j2 = It_Initial:N_DeltaTi:It_Final
% disp(['j2' ,' ', num2str(j2),' / ', num2str(It_Initial),' / ', num2str(It_Final)]);
% j2
nj2=nj2+1;AtwInteger=0.0;
for j1= 1:ntw
%% Emission Energy
aie= (HARMONICOMGA/OMGA1);
%% Gabor
xi = HARMONICOMGA*(a1(j1)-a1(j2));
W_xi = (Window^(-0.5)) * exp(-1.0*xi^2.0/(2.0*(Window)^2.00)) * exp(1i*xi)*(HARMONICOMGA^0.5);
W_tw = W_xi;
AtwInteger=AtwInteger+a2(j1)*W_tw*dti_AU;
end;
fprintf(fid1,'%32.26f %32.26f %32.26f\n',a1(j2)*OMGA1/(2.0*pi), aie, log10((abs(AtwInteger)^2.0)) );
fprintf(fid2,'%32.26f %32.26f %32.26f\n',a1(j2)*OMGA1/(2.0*pi), aie, ((abs(AtwInteger)^2.0)) );
end;
end; fclose(fid1); fclose(fid2);
%% ========================= 小窗时频分析 ========================= %%
%% ========================= 三维图:小窗时频分析(x~t,y~freq,z~intensity) ========================= %%
% surf(x,y,z)
load time_omga_yield.txt %读入时频分析数据
x=time_omga_yieldlog(:,1); %省略号换成你的x数据
y=time_omga_yieldlog(:,2);
z=time_omga_yieldlog(:,3);
'n_time_grids'
n_time_grids=floor((It_Final-It_Initial)/N_DeltaTi)
disp (['Number of time grids in wavelet-analysis : '])
nj2
figure(1003)
X=reshape(x,nj2,N_Shown_Energy_Range); % rephape数组 1维数组到 nt*n_freq 二维的数组 => mesh作图
Y=reshape(y,nj2,N_Shown_Energy_Range); % rephape数组 1维数组到 nt*n_freq 二维的数组
Z=reshape(z,nj2,N_Shown_Energy_Range); % rephape数组 1维数组到 nt*n_freq 二维的数组
mesh(X,Y,Z); % mesh作图
pcolor(X,Y,Z); % 类似于view001
%shading interp; % 伪彩色图
shading flat;
set(gca,'FontSize',12);
set(gca,'xtick',[ 0 1 2 3 4 5 6 7 8 9 10]); % x-grids
set(gca,'ytick',[5 7 9 11 13 15 17 19 21 23 25 27 29 31 ]); % y-grids
xlabel('time(o.c.)','FontSize',16), ylabel('harmonic order ','FontSize',16); % figure label
%colormap winter
colormap jet %winterrainbow
set(gca,'FontName','Helvetica');
set(gca,'Position',[.12 .17 .70 .74]);
set(gcf,'position',[100 100 780 460],'Units','centimeters');
%% ========================= 三维图:小窗时频分析(x~t,y~freq,z~intensity) ========================= %%
% mycmap = get(gcf,'Colormap')
% save('MyColormaps','mycmap')
评论0