不知道错在哪里了
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 = ;
% 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()-1;
P = 2^nextpow2(N);
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);
N=256;
fs=100;
f0=2.5;
a=0.7;
alpha=a*pi/2;
dt=1/fs;
t=(0:N-1)*dt;
df=fs/N;
f=(0:N-1)*df;
s=zeros(1,N)；
u=zeros(1,N);
i=sqrt(-1);
for m=1:N
s(m)=exp(-5*m*dt+i*cos(2*pi*f0*m*dt));
u(m)=sqrt((m*dt).^2+(m*df).^2).*sin(alpha+atan(m*dt./(m*df)));
end
F1=frft(s,a);
A1=F1'.*u;
B1=frft(A1,-a);
F2=frft(s,-a);
A2=F2'.*u;
B2=frft(A2,a);
dw=zeros(1,N);
w0=zeros(1,N);
w=zeros(1,N);
for m=1:N
dw(m)=(abs(B1').^2-abs(B2').^2)/(2*m*dt)/sin(2*alpha)/(abs(s).^2);
w(m)=dw(m)/(2*pi);
w0(m)=sin(2*pi*f0*m*dt);
end
plot(t,w,'r:',t,w0,'b-');
[ 本帖最后由 eight 于 2007-6-12 16:40 编辑 ] 原帖由 justwee 于 2007-6-12 16:25 发表 http://www.chinavib.com/forum/images/common/back.gif
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(nar ...
请给出出错信息，而不是要别人猜
回复 #2 eight 的帖子
我的程序是基于分数阶傅立叶变换的信号瞬时频率估计，运行显示如下错误：??? Undefined command/function 'sinc'.
Error in ==> frft>interp at 43
xint = fconv(y(1:2*N-1), sinc([-(2*N-3):(2*N-3)]'/2));
Error in ==> frft at 25
f = ;
Error in ==> QQ at 20
F1=frft(s,a); 原帖由 justwee 于 2007-6-12 19:01 发表 http://www.chinavib.com/forum/images/common/back.gif
我的程序是基于分数阶傅立叶变换的信号瞬时频率估计，运行显示如下错误：
??? Undefined command/function 'sinc'.
Error in ==> frft>interp at 43
xint = fconv(y(1:2*N-1), sinc([-(2*N-3):(2*N-3)]'/2 ...
sinc 是自定义的函数吗？要放在同一目录下
回复 #4 eight 的帖子
是sa函数，matlab 中应该有库函数的吧？ 原帖由 justwee 于 2007-6-12 20:15 发表 http://www.chinavib.com/forum/images/common/back.gif是sa函数，matlab 中应该有库函数的吧？
我这个版本自带 sinc 函数，不知道你的版本是否有，自己 help 一下吧
回复 #6 eight 的帖子
结果有了，但是和我的预想结果相差很多，为什么？回复 #6 eight 的帖子
红线应该和蓝线差不多才对的 原帖由 justwee 于 2007-6-12 20:23 发表 http://www.chinavib.com/forum/images/common/back.gif结果有了，但是和我的预想结果相差很多，为什么？
你的程序带有输入参数，旁人无法帮忙，自己检查一下吧
回复 #9 eight 的帖子
只有a可变，a是分数维傅立叶变换的维数。回复 #9 eight 的帖子
其它的不用考虑的。
页:
[1]