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')enda=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)));	endclose(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  endendclose(h);alpha=-N/2:N/2;f=-N/2:N/2;mesh(f,alpha,abs(S));xlabel('循環(huán)頻率alpha');ylabel('頻率f');zlabel('循環(huán)譜');