%两个非常接近的非线性变化的分数阶短时傅里叶变换和谱图分数阶谱图重排
clc;clear;close all;
fs=100;%采样频率
T = 16; %时宽
t =1/fs:1/fs:T; %所有采样点
%要处理的信号
% y = exp(2j*pi*(20*pi*t+5*pi*t.^2/2));
%等价于y = (cos(20*pi*t+5*pi*t.^2/2）*2*pi)
% % y = y/y(1);
% x=real(y');
y = cos(20*pi*t+2*pi*t.^2/2+2*pi*t.^3/40)+cos(24*pi*t+2*pi*t.^2/2+2*pi*t.^3/40);
y=y/y(1);
x=y';
%画出理论上的时频图
figure
finsh = (10+t+3*t.^2/40);
plot(finsh)
hold on
finsh = (12+t+3*t.^2/40);
plot(finsh)
title('理论计算得到的时频图');
%****特别重要：时频斜率归一化****
%通过归一化时频斜率求最佳阶数
%阶数是归一化时频斜率k的表达式角度 ? = arccot(-k)=arctan (1/k ) 然后计算阶次a = ? /（π/2）
%k = B/fs , B = mu*T mu是表达式中的斜率数值
% 原文链接：https://blog.csdn.net/weixin_42845306/article/details/120281437
dfinst=2*pi+(3*pi*t)/10; %每个点都对应一个mu，即表达式中的斜率数值
df_inst = dfinst/2/pi; %注意除去2pi
k = df_inst*T/fs %时频斜率归一化
alf= acot(-k)*2/pi %得到每个点的最佳阶数
%画出各个采样点最佳阶数的图像
figure
plot(alf)
title('各个采样点最佳阶数');
%初始化时间长度以及定义窗函数
[xrow,xcol] = size(x);
N =xrow;
bei = 4; %窗函数长度倍数
hname ='Blackman' ; %窗函数名字
hlength=floor(N/bei); %用矩形窗之后才能保持最佳阶数
hlength=hlength+1-rem(hlength,2); %强制转成奇数
t=1:xrow; h = tftb_window(hlength,hname); %用矩形窗之后才能保持最佳阶数
[hrow,hcol]=size(h); Lh=(hrow-1)/2; %Lh是窗函数长度（奇数）的一半
[trow,tcol] = size(t);
%经典的谱图和谱图重排
[sp,rsp,hat] = tfrrsp(x,1:N,N,h);
%求瑞丽熵
H_x_sp=MZJ_sang(sp)
H_x_rsp=MZJ_sang(rsp)
tfr= zeros(N,tcol); tf2= zeros(N,tcol); tf3= zeros (N,tcol);
indextfr= zeros(N,tcol);
%就是Th = t*h(t) Dh = h(t)对t求导就是调一下dwindow()函数
%先求出三个stft ：
% tfr： ht为窗函数
% tf2： Th = t*h(t)为窗函数
% tf3： Dh = h(t)对t求导为窗函数
Th=h.*[-Lh:Lh]'; Dh=dwindow(h);
for icol=1:tcol,
ti= t(icol);
tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau+Lh,N)+1;
norm_h=norm(h(Lh+1+tau));
tfr(indices,icol)=x(ti+tau,1).*conj( h(Lh+1+tau)) /norm_h;
tfr(:,icol)= MZJ_frft(tfr(:,icol),alf(icol)); %每一列做frft
indextfr(indices,icol) =x(ti+tau).*conj( h(Lh+1+tau)) /norm_h;
tf2(indices,icol)=x(ti+tau).*conj(Th(Lh+1+tau)) /norm_h;
tf3(indices,icol)=x(ti+tau).*conj(Dh(Lh+1+tau)) /norm_h;
num = 1;
indextfr(:,icol)= MZJ_frft(indextfr(:,icol),alf(icol));
tf2(:,icol)= MZJ_frft(tf2(:,icol),num); %每一列做frft
tf3(:,icol)= MZJ_frft(tf3(:,icol),num); %每一列做frft
end ;
%行变列,矩阵变成1列长度默认情况下是N*N
% %
f=(0.5*(0:N-1)/N)';
tfr=tfr(:);tf2=tf2(:);tf3=tf3(:); indextfr=indextfr(:);
%avoid_warn默认情况下就是1：N*N，因为没有一个数是0
avoid_warn=find(indextfr~=0);
%找到tfr矩阵中不等于0的下标
%例如：0012002401
%则 find(tfr~=0)得到
% 3 4 7 8 10
%为了坐标替换做准备
tf2(avoid_warn)=round(real(tf2(avoid_warn)./indextfr(avoid_warn)/1));
tf3(avoid_warn)=round(imag(N*tf3(avoid_warn)./indextfr(avoid_warn)/(2.0*pi)));
%谱图来了
tfr=abs(tfr).^2;
%把谱图三个东西都给重新规定大小
tfr=reshape(tfr,N,tcol);
tf2=reshape(tf2,N,tcol);
tf3=reshape(tf3,N,tcol);
%开始重排
rtfr= zeros(N,tcol);
Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-9*Ex;
for icol=1:tcol,
for jcol=1:N,
if abs(tfr(jcol,icol))>Threshold,
%t轴坐标替换real
icolhat= icol + tf2(jcol,icol);
icolhat=min(max(icolhat,1),tcol);
%f轴坐标替换imag
jcolhat= jcol - tf3(jcol,icol);
jcolhat=rem(rem(jcolhat-1,N)+N,N)+1;
%注意这时候不是0，
%因为这个jcolhat,icolhat前面生成的说不定就有值了
rtfr(jcolhat,icolhat)=rtfr(jcolhat,icolhat) + tfr(jcol,icol) ;
%这个累加本身才能来体现能量聚集
%刚开始的时候 rtfr(jcolhat,icolhat)是0
%这时候加的是 tfr(jcol,icol)加到(jcolhat,icolhat)位置
tf2(jcol,icol)=jcolhat + j * icolhat;
else
tf2(jcol,icol)=inf*(1+j);
rtfr(jcol,icol)=rtfr(jcol,icol) + tfr(jcol,icol) ;
end;
end;
end;
[tfrwv,t,f] = tfrwv(x);
H_x_tfrwv=MZJ_sang(tfrwv)
H_x_tfr=MZJ_sang(tfr)
H_x_rtfr=MZJ_sang(rtfr)
figure
imagesc(tfrwv);
xlabel('时间t');
ylabel('频率f') ;
title(['WVD']);
figure
imagesc(sp);
xlabel('时间t');
ylabel('频率f') ;
title(['窗函数',hname,'长度为N/',num2str(bei),'经典谱图']);
figure
imagesc(rsp);
xlabel('时间t');
ylabel('频率f');
title(['窗函数',hname,'长度为N/',num2str(bei),'经典谱图重排']);
% 谱图
figure
imagesc(t,f,tfr);
xlabel('时间t');
ylabel('频率f') ;
title(['窗函数',hname,'长度为N/',num2str(bei),'分数阶谱图']);
% 谱图重排
figure
imagesc(t,f,rtfr);
%mesh(t,f,rtfr);
xlabel('时间t');
ylabel('频率f');
title(['窗函数',hname,'长度为N/',num2str(bei),'分数阶谱图重排']);
figure
subplot(221)
mesh(t,f,sp);
xlabel('时间t');
ylabel('频率f') ;
title(['窗函数',hname,'长度为N/',num2str(bei),'经典谱图']);
subplot(222)
mesh(t,f,rsp);
xlabel('时间t');
ylabel('频率f') ;
title(['窗函数',hname,'长度为N/',num2str(bei),'经典谱图重排']);
subplot(223)
mesh(t,f,tfr);
xlabel('时间t');
ylabel('频率f') ;
title(['窗函数',hname,'长度为N/',num2str(bei),'分数阶谱图']);
subplot(224)
mesh(t,f,rtfr);
xlabel('时间t');
ylabel('频率f') ;
title(['窗函数',hname,'长度为N/',num2str(bei),'分数阶谱图重排']);
