傅立叶变换是将观看角度从时域转变到频域，分数阶傅立叶变换就是以观看时频面的角度去旋转时频面的坐标轴，然后再从观察频域的角度去分析信息。
分数阶傅立叶变换多出来的一个算子就是这个旋转的角度。这个旋转角度以分数的形式呈现，取值是0-1，当取1时就等同于傅立叶变换。
将信息进行分数阶傅立叶变换的原因在于：大部分信息都是非平稳信号，仅仅用傅立叶变换不足以分析其显著特征，运用分数阶傅立叶变换主要是能选取信息最集中的角度去分析，也就是在不同的分数阶得到的结果中选取幅值最大的那个结果，那么这个结果所存在的那个分数阶就是最优阶次。
FRFT定义
一般地，信号
x
(
t
)
x(t)
x(t)，其
p
p
p阶分数阶 Fourier变换可表示为
F
p
x
F^px
Fpx 或
X
ϕ
(
u
)
X_\phi(u)
Xϕ(u)，前者可以解释为算子
F
p
F^p
Fp作用于函数
x
(
t
)
x(t)
x(t)，后者可以解释为分数阶 Fourier 变换的结果在
u
u
u 域上，
ϕ
=
p
π
/
2
\phi = p\pi/2
ϕ=pπ/2 为阶数
p
p
p 对应的旋转角度，即新变量
u
−
u-
u−轴相对于时间
t
−
t-
t−轴逆时针转过的角度。具体地，FRFT 定义为:
F
p
x
=
X
ϕ
(
u
)
=
∫
−
∞
∞
x
(
t
)
K
ϕ
(
t
,
u
)
d
t
(1)
F^px=X_\phi(u)=\int_{-\infty}^{\infty} x(t)K_\phi(t,u) dt \tag{1}
Fpx=Xϕ(u)=∫−∞∞x(t)Kϕ(t,u)dt(1)
其中
K
ϕ
(
t
,
u
)
K_\phi(t,u)
Kϕ(t,u) 为分数阶 Fourier 变换的核函数，定义如下：
K
ϕ
(
t
,
u
)
=
{
A
ϕ
e
x
p
(
j
t
2
+
u
2
2
c
o
t
ϕ
−
j
u
t
c
s
c
ϕ
)
,
f
o
r
ϕ
≠
n
π
δ
(
t
−
u
)
,
f
o
r
ϕ
=
2
n
π
δ
(
t
+
u
)
,
f
o
r
ϕ
=
(
2
n
±
1
)
π
(2)
K_\phi(t,u)=\begin{cases} A_\phi exp(j^{\frac{t^2+u^2}{2}}cot\phi-j\space ut\space csc\phi) ,for\phi\neq n\pi\\ \delta(t-u),for\phi=2n\pi\\\delta(t+u),for\phi=(2n\pm1)\pi \end{cases}\tag{2}
Kϕ(t,u)=⎩⎪⎨⎪⎧Aϕexp(j2t2+u2cotϕ−j ut cscϕ),forϕ=nπδ(t−u),forϕ=2nπδ(t+u),forϕ=(2n±1)π(2)其中
A
ϕ
=
1
−
j
c
o
t
ϕ
2
π
A_\phi=\sqrt{\frac{1-j\space cot\phi}{2\pi}}
Aϕ=2π1−j cotϕ为幅度因子。
FRFT 的核函数和逆变换
分数阶 Fourier 变换的核函数具有如下性质：
K
ϕ
(
t
,
u
)
=
K
ϕ
(
u
,
t
)
(3)
K_\phi(t,u)=K_\phi(u,t)\tag{3}
Kϕ(t,u)=Kϕ(u,t)(3)
K
−
ϕ
(
t
,
u
)
=
K
ϕ
∗
(
t
,
u
)
(4)
K_{-\phi}(t,u)=K^*_\phi(t,u)\tag{4}
K−ϕ(t,u)=Kϕ∗(t,u)(4)
K
ϕ
(
−
t
,
u
)
=
K
ϕ
(
t
,
−
u
)
(5)
K_\phi(-t,u)=K_\phi(t,-u)\tag{5}
Kϕ(−t,u)=Kϕ(t,−u)(5)
∫
−
∞
∞
K
ϕ
(
t
,
u
)
K
φ
(
u
,
z
)
d
u
=
K
ϕ
+
φ
(
t
,
z
)
(6)
\int_{-\infty}^{\infty} K_\phi(t,u)K_\varphi(u,z)du=K_{\phi+\varphi}(t,z)\tag{6}
∫−∞∞Kϕ(t,u)Kφ(u,z)du=Kϕ+φ(t,z)(6)
∫
−
∞
∞
K
ϕ
(
t
,
u
)
K
ϕ
∗
(
t
,
u
′
)
d
t
=
δ
(
u
−
u
′
)
(7)
\int_{-\infty}^{\infty} K_\phi(t,u)K^*_\phi(t,u^\prime)dt=\delta(u-u^\prime)\tag{7}
∫−∞∞Kϕ(t,u)Kϕ∗(t,u′)dt=δ(u−u′)(7)应用这些性质，容易得到结论：旋转角度为
ϕ
\phi
ϕ的分数阶 Fourier变换的逆变换即为旋转角度为
−
ϕ
-\phi
−ϕ 的分数阶 Fourier 变换。对于函数
x
(
t
)
x(t)
x(t)，用公式可表示为
x
(
t
)
=
∫
−
∞
∞
X
ϕ
(
u
)
K
−
ϕ
(
t
,
u
)
d
u
(8)
x(t)=\int_{-\infty}^{\infty}X_\phi(u)K_{-\phi}(t,u)du \tag{8}
x(t)=∫−∞∞Xϕ(u)K−ϕ(t,u)du(8)显然，当阶数 p = 1 时，有
ϕ
=
π
/
2
\phi= \pi/2
ϕ=π/2，式 (2) 即为传统的 Fourier 变换核，式 (8)也即为传统的 Fourier 逆变换。因此，分数阶 Fourier 变换是一种广义的 Fourier 变换。由核函数的定义及其性质，分数阶 Fourier 变换是一种线性变换，其基底是一组正交的 chirp 函数，即线性调频的复指数函数。对于不同的
u
u
u值，核函数的区别只有一个与
u
u
u有关的时延和相位：
K
ϕ
(
t
,
u
)
=
e
−
j
u
2
2
t
a
n
ϕ
K
ϕ
(
t
−
u
s
e
c
ϕ
,
0
)
(9)
K_\phi(t,u)=e^{-j\frac{u^2}{2}tan\phi}K_\phi(t-u\space sec \space \phi,0)\tag{9}
Kϕ(t,u)=e−j2u2tanϕKϕ(t−u sec ϕ,0)(9)
FRFT 的性质
由其核函数，我们可以得到分数阶 Fourier 变换的许多性质：
1.线性变换特性：
F
p
[
∑
n
c
n
x
n
(
t
)
]
=
∑
n
c
n
[
F
p
x
n
(
t
)
]
(10)
F^p[\begin{matrix} \sum_{n}c_nx_n(t) \end{matrix} ]=\begin{matrix} \sum_{n}c_n[F^p x_n(t)] \end{matrix}\tag{10}
Fp[∑ncnxn(t)]=∑ncn[Fpxn(t)](10)2.旋转相加性：
F
p
1
F
p
2
=
F
p
1
+
p
2
(11)
F^{p1}F^{p2}=F^{p1+p2}\tag{11}
Fp1Fp2=Fp1+p2(11)3.阶数为整数 n 的分数阶 Fourier 变换相当于传统 Fourier 变换的 n 次幂，即
F
n
=
(
F
)
n
(12)
F^{n}=(F)^{n}\tag{12}
Fn=(F)n(12)由于传统 Fourier 变换具有周期性，即
F
4
=
F
0
F4 = F0
F4=F0，因此还有
F
n
=
F
m
,
f
o
r
n
=
m
o
d
(
m
,
4
)
(13)
F^n=F^m,for\space n=mod(m,4)\tag{13}
Fn=Fm,for n=mod(m,4)(13)4.时移特性：
{
F
p
[
x
(
t
−
τ
)
]
}
(
u
)
=
e
j
π
τ
2
s
i
n
ϕ
c
o
s
ϕ
−
j
2
π
u
τ
s
i
n
ϕ
X
ϕ
(
u
−
τ
c
o
s
ϕ
)
(14)
\{F^p[x(t-\tau)]\}(u)=e^{j\pi\tau^2sin\phi cos\phi-j2\pi u\tau\space sin\phi}X_\phi(u-\tau cos\phi)\tag{14}
{Fp[x(t−τ)]}(u)=ejπτ2sinϕcosϕ−j2πuτ sinϕXϕ(u−τcosϕ)(14)5.频移特性：
{
F
p
[
e
j
2
π
ξ
t
x
(
t
)
]
}
(
u
)
=
e
−
j
π
ξ
2
s
i
n
ϕ
c
o
s
ϕ
−
j
2
π
u
ξ
c
o
s
ϕ
X
ϕ
(
u
−
ξ
s
i
n
ϕ
)
(15)
\{F^p[e^{j2\pi\xi t}x(t)]\}(u)=e^{-j\pi\xi^2sin\phi cos\phi-j2\pi u\xi\space cos\phi}X_\phi(u-\xi sin\phi)\tag{15}
{Fp[ej2πξtx(t)]}(u)=e−jπξ2sinϕcosϕ−j2πuξ cosϕXϕ(u−ξsinϕ)(15)6.Parseval 准则：
∫
−
∞
∞
f
∗
(
u
)
g
(
u
)
d
u
=
∫
−
∞
∞
f
ϕ
∗
(
u
)
g
ϕ
(
u
)
d
u
(16)
\int_{-\infty}^{\infty}f^*(u)g(u)du=\int_{-\infty}^{\infty}f_\phi^*(u)g_\phi(u)du \tag{16}
∫−∞∞f∗(u)g(u)du=∫−∞∞fϕ∗(u)gϕ(u)du(16)
LFM信号的FRFT
由式 (8) 可知，对任意信号 x ( t ) x(t) x(t)，其分数阶 Fourier 变换 X ϕ ( u ) X\phi(u) Xϕ(u) 可以看作为 x ( t ) x(t) x(t) 在以变换核 K − ϕ ( u , t ) K−\phi(u, t) K−ϕ(u,t)为基底的函数空间的展开。由于分数阶 Fourier 变换核是 u u u域上的一组正交 chirp 基，因此 LFM 信号在一个合适的分数阶 Fourier 域中能量可以被汇聚为一个冲激函数。具体地，考虑如下有限长度的 LFM 信号： s ( t ) = r e c t ( t ) e x p ( j ( 2 π f 0 t + ϕ k t 2 + φ 0 ) ) (17) s(t)=rect(t)exp(j(2\pi f_0t+\phi kt^2+\varphi_0))\tag{17} s(t)=rect(t)exp(j(2πf0t+ϕkt2+φ0))(17)其中 f 0 , k , φ 0 f_0,k,\varphi_0 f0,k,φ0分别为起始频率、调频斜率和初始相位，终止频率为 f 1 = f 0 + k T f1= f0+ kT f1=f0+kT ， T T T为信号持续时间，带宽即为 B = f 1 − f 0 = k T B = f_1− f_0= kT B=f1−f0=kT；当 t ∈ [ 0 , T ] t\in[0, T ] t∈[0,T]时，归一化的矩形窗函数 r e c t ( t ) = 1 / T rect (t) = 1/T rect(t)=1/T，此外值为 0 0 0。瞬时频率，即相位对时间的导数： f L ( t ) = 1 2 π d d t ( 2 π f 0 t + π k t 2 + φ 0 ) = f 0 + k t (18) f_L(t)=\frac{1}{2\pi}\frac{d}{dt}(2\pi f_0t+\pi kt^2+\varphi_0)=f_0+kt\tag{18} fL(t)=2π1dtd(2πf0t+πkt2+φ0)=f0+kt(18)在时频平面上是关于t的线性函数，如下图。
最佳旋转角度与 LFM 信号调频斜率的关系是
ϕ
⋆
=
−
a
r
c
t
a
n
(
1
k
)
(19)
\phi^⋆=-arctan(\frac{1}{k})\tag{19}
ϕ⋆=−arctan(k1)(19)即
k
=
−
c
o
t
ϕ
⋆
(20)
k=-cot\phi^⋆\tag{20}
k=−cotϕ⋆(20)
参考代码为：论文Adhemar Bultheel and H´ector E. Mart´ınez Sulbaran.frft:Computation of the Fractional Fourier Transform,2004.FRFT代码链接
function Faf = frft(f, a)
% The fast Fractional Fourier Transform
% input: f = samples of the signal
% a = fractional power
% output: Faf = fast Fractional Fourier transform
error(nargchk(2, 2, nargin));
f = f(:);
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;
sN = sqrt(N);
a = mod(a,4);
% do special cases
if (a==0), Faf = f; return; end;
if (a==2), Faf = flipud(f); return; end;
if (a==1), Faf(shft,1) = fft(f(shft))/sN; return; end
if (a==3), Faf(shft,1) = ifft(f(shft))*sN; return; end
% reduce to interval 0.5 < a < 1.5
if (a>2.0), a = a-2; f = flipud(f); end
if (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end
if (a<0.5), a = a+1; f(shft,1) = ifft(f(shft))*sN; end
% the general case for 0.5 < a < 1.5
alpha = a*pi/2;
tana2 = tan(alpha/2);
sina = sin(alpha);
f = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)];
% chirp premultiplication
chrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);
f = chrp.*f;
% chirp convolution
c = pi/N/sina/4;
Faf = fconv(exp(i*c*(-(4*N-4):4*N-4)'.^2),f);
Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);
% chirp post multiplication
Faf = chrp.*Faf;
% normalizing constant
Faf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);
%%%%%%%%%%%%%%%%%%%%%%%%%
function xint=interp(x)
% sinc interpolation
N = length(x);
y = zeros(2*N-1,1);
y(1:2:2*N-1) = x;
xint = fconv(y(1:2*N-1), sinc([-(2*N-3):(2*N-3)]'/2));
xint = xint(2*N-2:end-2*N+3);
%%%%%%%%%%%%%%%%%%%%%%%%%
function z = fconv(x,y)
% convolution by fft
N = length([x(:);y(:)])-1;
P = 2^nextpow2(N);
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);
同时贴出Python版本的主函数代码。
"""
Module to calculate the fast fractional fourier transform.
"""
from __future__ import division
import numpy
import scipy
import scipy.signal
def frft(f, a):
"""
Calculate the fast fractional fourier transform.
Parameters
----------
f : numpy array
The signal to be transformed.
a : float
fractional power
Returns
-------
data : numpy array
The transformed signal.
References
---------
.. [1] This algorithm implements `frft.m` from
https://nalag.cs.kuleuven.be/research/software/FRFT/
"""
ret = numpy.zeros_like(f, dtype=numpy.complex)
f = f.copy().astype(numpy.complex)
N = len(f)
shft = numpy.fmod(numpy.arange(N) + numpy.fix(N / 2), N).astype(int)
sN = numpy.sqrt(N)
a = numpy.remainder(a, 4.0)
# Special cases
if a == 0.0:
return f
if a == 2.0:
return numpy.flipud(f)
if a == 1.0:
ret[shft] = numpy.fft.fft(f[shft]) / sN
return ret
if a == 3.0:
ret[shft] = numpy.fft.ifft(f[shft]) * sN
return ret
# reduce to interval 0.5 < a < 1.5
if a > 2.0:
a = a - 2.0
f = numpy.flipud(f)
if a > 1.5:
a = a - 1
f[shft] = numpy.fft.fft(f[shft]) / sN
if a < 0.5:
a = a + 1
f[shft] = numpy.fft.ifft(f[shft]) * sN
# the general case for 0.5 < a < 1.5
alpha = a * numpy.pi / 2
tana2 = numpy.tan(alpha / 2)
sina = numpy.sin(alpha)
f = numpy.hstack((numpy.zeros(N - 1), sincinterp(f), numpy.zeros(N - 1))).T
# chirp premultiplication
chrp = numpy.exp(-1j * numpy.pi / N * tana2 / 4 *
numpy.arange(-2 * N + 2, 2 * N - 1).T ** 2)
f = chrp * f
# chirp convolution
c = numpy.pi / N / sina / 4
ret = scipy.signal.fftconvolve(
numpy.exp(1j * c * numpy.arange(-(4 * N - 4), 4 * N - 3).T ** 2),
f
)
ret = ret[4 * N - 4:8 * N - 7] * numpy.sqrt(c / numpy.pi)
# chirp post multiplication
ret = chrp * ret
# normalizing constant
ret = numpy.exp(-1j * (1 - a) * numpy.pi / 4) * ret[N - 1:-N + 1:2]
return ret
def ifrft(f, a):
"""
Calculate the inverse fast fractional fourier transform.
Parameters
----------
f : numpy array
The signal to be transformed.
a : float
fractional power
Returns
-------
data : numpy array
The transformed signal.
"""
return frft(f, -a)
def sincinterp(x):
N = len(x)
y = numpy.zeros(2 * N - 1, dtype=x.dtype)
y[:2 * N:2] = x
xint = scipy.signal.fftconvolve(
y[:2 * N],
numpy.sinc(numpy.arange(-(2 * N - 3), (2 * N - 2)).T / 2),
)
return xint[2 * N - 3: -2 * N + 3]
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
