马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我想把这段程序换成C语言,但由于没学过matlab,所有有些地方搞不明白,希望matlab高手帮忙说明和注释一下,万分感谢,在线等待!
function S=fft_accumulation_method(x,y,N,a,g,L)
%
% FFT_ACCUMULATION_METHOD
% estimate the spectral correlation density using the FFT
% accumulation method
% Reference:
% Roberts R. S., Brown, W. A. and Loomis, H. H.
% "Computationally efficient algorithms for cyclic
% spectral analysis" IEEE Signal Processing Magazine
% 8(2) pp38-49 April 1991
%
% Input parameters are:
% x,y signals
% N length of time window used for estimating frequency
% segments (should be power of 2)
% a window used for smoothing segments
% g window for smoothing correlation
% L decimation factor
%
% USAGE S=fft_accumulation_method(x,y,N,a,g,L)
%
% e.g s=fft_accumulation_method(s1,s1,64,'hamming','hamming',1)
if ~exist('L')
L=1;
end
lx=length(x);
if (length(y)~=lx)
error('Time series x and y must be same length')
end
n=0:floor((lx-N)/L);
ln=length(n);
a=feval(a,N)';
g=feval(g,ln)';
g=g/sum(g);
a=a/sum(a);
X=zeros(N,ln);
Y=zeros(N,ln);
Ts=1/N;
for i=1:ln
n_r=n(i)*L+(1:N);
X(:,i)=fftshift(fft(a.*x(n_r)))';
Y(:,i)=fftshift(conj(fft(a.*y(n_r))))';
end
lnd2=floor(ln/2);
lnd4=floor(ln/4)+1;
ln3d4=lnd4+lnd2-2;
S=zeros(2*N-1,lnd2*(N+1));
for alpha=-N/2+1:N/2-1
for f=-N/2:N/2-1
f1=f+alpha;
f2=f-alpha;
if ((abs(f1)<N/2)&(abs(f2)<N/2))
f1=f1+N/2;
f2=f2+N/2;
fsh=fftshift(fft(g.*X(f1,:).*Y(f2,:)));
S(2*f+N,(alpha+N/2)*lnd2+(1:lnd2-1))=fsh(1,lnd4:ln3d4);
end
f1=f+alpha;
f2=f-alpha+1;
if ((abs(f1)<N/2)&(abs(f2)<N/2))
f1=f1+N/2;
f2=f2+N/2;
fsh=fftshift(fft(g.*X(f1,:).*Y(f2,:)));
S(2*f+N+1,(alpha+N/2)*lnd2+(1:lnd2-1)-lnd4+1)=fsh(1,lnd4:ln3d4);
end
end
end |