fft-example.zip_快速傅里叶变换例子

% 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')