function S=cyclic_cross_periodogram(x,y,N,M,a,g)
%close all;
%clear all;
% CYCLIC_CROSS_PERIODOGRAM
%
% calculate the spectral correlation density function
% of signals x and y using the cyclic cross periodogram
% algorithm
%
% Input parameters are:
% x,y signals
% N length of time window used for estimating frequency
% segments
% a window used for smoothing segments
% g window for smoothing correlation
%
% USAGE S=cyclic_cross_periodogram(x,y,N,a,g)
%
% e.g s=cyclic_cross_periodogram(s1,s1,64,'hamming','hamming')
lx=length(x);
if (length(y)~=lx)
error('Time series x and y must be same length')
end
a=feval(a,N)';
g=feval(g,M)';
g=g/sum(g);
a=a/sum(a);
X=zeros(2*N+1);
Y=zeros(2*N+1);
Ts=1/N;
a=waitbar(0,'f iteration');
for f=-N:N
waitbar((f+N+1)/(2*N+1))
xf=x.*exp(-j*2*pi*f*(0:lx-1)*Ts);
yf=y.*exp(-j*2*pi*f*(0:lx-1)*Ts);
n_r=N:-1:1;
X(f+N+1)=sum(a.*xf(n_r));
Y(f+N+1)=conj(sum(a.*yf(n_r)));
end
close(a)
h=waitbar(0,'alpha iteration');
S=zeros(N+1,N+1);
for alpha=-floor(N/2):floor(N/2)
waitbar((alpha+N/2)/N);
for f=-floor(N/2):floor(N/2)
f1=f+floor(alpha/2)+(floor(-((M-1)/2)):floor((M-1)/2));
f2=f-floor(alpha/2)+(floor(-((M-1)/2)):floor((M-1)/2));
if ((abs(f1)<N)&(abs(f2)<N))
S(f+floor(N/2)+1,floor(N/2)+alpha+1)=sum(g.*X(f1+N+1).*Y(f2+N+1))/M;
end
end
end
close(h);
alpha=-N/2:N/2;
f=-N/2:N/2;
mesh(f,alpha,abs(S));
xlabel('循环频率alpha');
ylabel('频率f');
zlabel('循环谱');
cyclic_cross_periodogram.rar_cyclostationary_三维循环谱_三维谱_循环平稳谱_谱相关
版权申诉
5星 · 超过95%的资源 25 浏览量
2022-07-14
07:03:03
上传
评论
收藏 1KB RAR 举报
JonSco
- 粉丝: 72
- 资源: 1万+