function Faf = frft(f, a)
narginchk(2, 2); %用来在M文件中检查你的输入的变量数是否合格，你可以设置第一个参数low，即最小的传递参数的数目，第二个参数high，最大的传递参数个数。这里两个都是2，则必须是传递两个参数。如果参数个数满足大于low，小于high，则返回一个空矩阵，否则返回错误消息，由error（msg）显示出来。
f = f(:); %将矩阵f的每一列元素堆积起来，成为一个列向量
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1; %总体是将右边的一半数移到左边；fix是舍零取整数的意思，就是1.9或1.1也是1；
sN = sqrt(N);
a = mod(a,4); %阶数是一个周期为4的函数；mod函数和rem函数的区别。在这里用一个方法判断，取余遵循尽可能让余数的绝对值小的原则，取模遵循尽可能让商小的原则，
% 特殊情况
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;且frft是关于阶数的周期函数
if (a>2.0), a = a-2; f = flipud(f); end %表示将参与frft的信号变为原始信号的上下反转（alpha = a*pi/2）
if (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end %即将参与frft的信号变为原始信号f的的fft
if (a<0.5), a = a+1; f(shft,1) = ifft(f(shft))*sN; end %即将参与的frft信号变为原始信号f的ifft
% 一般情况 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)]; %前后各插入N个0，将f长度变为4N
% 第一步：原信号与一个线性调频信号chrp的左乘
chrp = exp(-1i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);
f = chrp.*f;
% 第二步：左乘后的结果与另一个函数卷积
c = pi/N/sina/4;
Faf = fconv(exp(1i*c*(-(4*N-4):4*N-4)'.^2),f);
Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi); %去卷积结果的中间1/3总长度的数据点数
% 第三步：卷积结果再与chrp相乘
Faf = chrp.*Faf;
% 归一化常数
Faf = exp(-1i*(1-a)*pi/4)*Faf(N:2:end-N+1);
function xint=interp(x) %基于FFT的内插法
% 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); %取卷积结果的中间1/3总长度的数据,end表示最后一行
function z = fconv(x,y) % 圆周卷积
% convolution by fft
N = length([x(:);y(:)])-1;
P = 2^nextpow2(N); %取最接近数据长度的2的整数次方，比如N=1021，执行n=2^nextpow2(N)后，n=2^10；可以使fft更快
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);
评论0