各种傅里叶变换

clc,clear,close all;
N=1024;
sig1=fmlin(N,0.15,0.25);
%%暂态信号2
sig2=fmlin(N,0.1,0.2);
sig=sig1+sig2;
%这个是估算瞬时频率
figure(3);
[f,t]=instfreq(sig1);
plot(t,f);
%短时傅里叶变换
%就是这个tfrstft
stft=MZJ_stft( sig );
figure(1);
imagesc( abs(stft) ); 
figure(2);
mesh( abs(stft) ); 

clc;clear
%知识点：tfr源代码就是stft的最后做完fft之后平方一次

% fs=300;
% t =1/fs:1/fs:15;
N = 128;
sig1 = fmlin(N,0,0.3);
sig2 = fmlin(N,0.2,0.4);
sig = sig1+sig2;
% 长度为63的高斯窗时,计算得到的谱图
% 由于两个线性调频信号相隔较远,所以取长的频域平滑窗,可以消除谱图中的交叉项。
%设置窗函数
h1 = tftb_window(23,'Gauss');
% h1=gausswin(23);
[sp,t,f] = tfrsp(sig,1:N,N,h1);
%为了保持f是递增的就加上fftshift（）
contour(t,fftshift(f),sp); 
xlabel('时间t');
ylabel('频率f') ;
title('sp'); 
figure(2);
%这个窗函数的长度只能是pow（2，n）-1
h2=gausswin(63);
[sp,t,f] = tfrsp(sig,1:N,N,h2);
% imagesc(t,f(1:length(f)/2),sp);
contour(t,fftshift(f),(sp)); 
xlabel('时间t');
ylabel('频率f') ;
title('sp'); 

1.窗函数tftb_window

function h=tftb_window(N,name,param,param2);
%tftb_window	Window generation.
%	N      : length of the window
%	NAME   : name of the window shape (default : Hamming)
%	PARAM  : optional parameter
%	PARAM2 : second optional parameters
%
%	Possible names are :
%	'Hamming', 'Hanning', 'Nuttall',  'Papoulis', 'Harris',
%	'Rect',    'Triang',  'Bartlett', 'BartHann', 'Blackman'
%	'Gauss',   'Parzen',  'Kaiser',   'Dolph',    'Hanna'.
%	'Nutbess', 'spline',  'Flattop'
%       h=tftb_window(256,'Gauss',0.005); 

%%
%探究wvd的实现代码
%问题就是为什么  indices就是要做fft的函数对应的时间轴
%取到了1后面一部分，N后面一部分还是对称的
clc;clear
N = 128;
sig1 = 2*fmlin(N,0,0.3);
sig2 =3*fmlin(N,0.2,0.4);
sig = sig1+sig2;

x = fmlin(N,0,0.3);

% xcol必须为一
%xrow就是N
[xrow,xcol] = size(x);
 t=1:xrow; N=xrow ;
[trow,tcol] = size(t);

tfr= zeros (N,tcol);  
for icol=1:tcol,
    %这个语句就是说 从1到tcol默认为N
%  ti= t(icol);
 ti= icol;
 %round用于舍入到最接近的整数。
 %找到的taumax是距离0或者N最小的值
 %例如ti = 4
 taumax=min([ti-1,xrow-ti,round(N/2)-1]);
 %taumax=3
 tau=-taumax:taumax; 
 %-3：3
  %r = rem(a,b) 返回 a 除以 b 后的余数
  %indices就是要做fft的函数对应的时间轴
 indices= rem(N+tau,N)+1;
 %取到了【126,127,128,1,2,3,4】
 
 %以ti为中心向两边扩散，直到碰到0，两边长度相同
 tfr(indices,icol) = x(ti+tau,1) .* conj(x(ti-tau,xcol));
 %这个就是求那个做fft的函数
 %126 ：x(1)*conj(x(7))
 %127：x(2)*conj(x(6))
 %128：x(3)*conj(x(5))
  %1：x(4)*conj(x(4))
  %开始对称
  %2：x(5)*conj(x(3))
  

 tau=round(N/2); 
 if (ti<=xrow-tau)&(ti>=tau+1),
  tfr(tau+1,icol) = 0.5 * (x(ti+tau,1) * conj(x(ti-tau,xcol))  + ...
                           x(ti-tau,1) * conj(x(ti+tau,xcol))) ;
 end;
end; 
%就是对每一列做fft，而indices就是要做fft的函数对应的时间轴
% tfr= fft(tfr(1:round(N/2),1:N)); 
tfr= fft(tfr); 
% % f1 = tfr(:,8)
% % f2 = fft(f1);
% if (xcol==1), tfr=real(tfr); end ;
f=(0.5*(0:N-1)/N)';
contour(tfr);

1.stft实现代码

%%
clc;clear
N = 64;
sig1 = 2*fmlin(N,0,0.3);
sig2 =3*fmlin(N,0.2,0.4);
sig = sig1+sig2;
x = fmlin(N,0,0.3);
% xcol必须为一
%xrow就是N
[xrow,xcol] = size(x);
 N=xrow ;
hlength=floor(N/4);
%强制转成奇数
hlength=hlength+1-rem(hlength,2);
 t=1:xrow; h = tftb_window(hlength); 


[trow,tcol] = size(t);


[hrow,hcol]=size(h);
%hcol也是长度
%Lh是窗函数长度（奇数）的一半
Lh=(hrow-1)/2; 
%N/2

h=h/norm(h);
plot(h)
%窗函数幅度归一化
tfr= zeros (N,tcol);  
for icol=1:tcol,%每行对应的一个点+[-N/2,N/2]上做fft
 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;    
 %把范围限制在[1,N]rem(N+tau+Lh,N)+1,把全部的顺序排列
 %rem(N+tau,N)+1,把一半在最左边，一半在最右边
indices
tfr(indices,icol)=x(ti+tau,1).*conj(h(Lh+1+tau));
end
%就是对每一列做fft，而indices就是要做fft的函数对应的时间轴
% tfr= fft(tfr(1:round(N/2),1:N)); 
tfr= fft(tfr); 
% % f1 = tfr(:,8)
% % f2 = fft(f1);
% if (xcol==1), tfr=real(tfr); end ;
f=(0.5*(0:N-1)/N)';
contour(t,f,abs(tfr)); 

3. rsp

​                               

  

下面公式是错的但是最后一行比较重要
 

%%
%谱图重排
clc;clear
fs=300;
t =1/fs:1/fs:15;
y = sin(20*pi*t+5*pi*t.^2/2+2*pi*t.^3/9-pi*t.^4/80)+sin(24*pi*t+5*pi*t.^2/2+2*pi*t.^3/9-pi*t.^4/80);

N = 1280;
sig1 = fmlin(N,0,0.3);
sig2 = fmlin(N,0.2,0.4);
sig = sig1+sig2;

%谱图
figure(1)
[sp,rsp,hat]=tfrrsp(sig);
imagesc(sp); 
xlabel('时间t');
ylabel('频率f') ;
title('谱图'); 

%谱图重排
figure(2)
imagesc(rsp);
xlabel('时间t');
ylabel('频率f');
title('谱图重排'); 
%%
%谱图重排底层原理
clc;clear

N = 128;
sig1 = fmlin(N,0,0.3);
sig2 = fmlin(N,0.2,0.4);
x = sig1+sig2;
[xrow,xcol] = size(x);

hlength=floor(N/4);
hlength=hlength+1-rem(hlength,2);
 t=1:xrow; h = tftb_window(hlength);
 [hrow,hcol]=size(h); Lh=(hrow-1)/2; 
 [trow,tcol] = size(t);


 
tfr= zeros(N,tcol); tf2= zeros(N,tcol); tf3= 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,N)+1;
 norm_h=norm(h(Lh+1+tau));
 tfr(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;
end ;
tfr=fft(tfr);
tf2=fft(tf2);
tf3=fft(tf3);
%得到了3个stft
%行变列,矩阵变成1列长度默认情况下是N*N
tfr=tfr(:);tf2=tf2(:);tf3=tf3(:);
%avoid_warn默认情况下就是1：N*N，因为没有一个数是0
avoid_warn=find(tfr~=0);
%找到tfr矩阵中不等于0的下标
%例如：0012002401
%则 find(tfr~=0)得到
%  3  4  7  8  10
%为了坐标替换做准备
tf2(avoid_warn)=round(real(  tf2(avoid_warn)./tfr(avoid_warn)/1));
tf3(avoid_warn)=round(imag(N*tf3(avoid_warn)./tfr(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-6*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;
%谱图
figure(1)
imagesc(tfr); 
xlabel('时间t');
ylabel('频率f') ;
title('谱图'); 
%谱图重排
figure(2)
imagesc(rtfr);
xlabel('时间t');
ylabel('频率f');
title('谱图重排');