matlab中分数阶线性系统的bode图程序看不懂，求高手指点，越详细越好，多谢啦！

function H=f_bode(a,eta,b,gamma,w,Gc)

%分数笑旦阶线性系统的bode图程序

%a,b,eta,gamma均为整数，w实值矢量，Gc权重，输出H是frd模型，例碰腊扰如：

%a=3b=1eta=3gamma=2freq = logspace(1,2)w = 0.05*(freq)f_bode(a,eta,b,gamma,w)

%a=3b=1eta=3gamma=2freq = linspace(0.01,0.2)w = 0.05*(freq)f_bode(a,eta,b,gamma,w)

H=zeros(size(w))

for i=1:length(w)

P=b*((sqrt(-1)*w(i)).^gamma.')

Q=a*((sqrt(-1)*w(i)).^eta.')

H(i)=P/Q

end

H=frd(H,w)%生成频率响应局野数据模型

if nargin==6,

H=H*Gc

end

if nargout==0,

bode(H)

end

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 = freturnend

if (a==2), Faf = flipud(f)returnend

if (a==1), Faf(shft,1) = fft(f(shft))/sNreturnend

if (a==3), Faf(shft,1) = ifft(f(shft))*sNreturnend

% reduce to interval 0.5 <a <1.5

if (a>2.0), a = a-2f = flipud(f)end

if (a>1.5), a = a-1f(shft,1) = fft(f(shft))/sNend

if (a<0.5), a = a+1f(shft,1) = ifft(f(shft))*sNend

% 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)

用MATLAB 实现傅里叶变换: 用户任意输入一个函数，然后，输出函袭搏历数的傅里叶变换函数，然后输出振幅频率 。拍搜 x=sin(2*pi*t)%任意输入一个函数。 y=fft(x)%傅里银纤叶变换函数。 plot(abs(y))%振幅频率。

---