www.pudn.com > cyclic2.rar > cyclic2.m, change:2005-12-28,size:2744b


clear 
clc 
% bark1=[1,1,1,1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,1,-1,-1,1,1,1,-1,1,-1,1,-1,1,1,1,-1,1,-1,1,-1,1,1,1,-1,1,-1,1]; 
bark1=[-1,1,-1,1,-1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,-1,1,-1,1,-1,1,-1,-1,1,-1,-1,1,-1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,1,1,-1,1,-1,1,1,1,1,1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,1,-1,-1,1,1,1,-1,1,-1,1,-1,1,1,1,-1,1,-1,1,-1,1,1,1,-1,1,-1,1]; 
bark=[bark1 bark1 bark1]; 
% barklen=length(bark) 
 
fs=1;%采样频率 
Ts=1/fs; 
fb=1/8; 
Tb=1/fb; 
% N=13*fs/fb; 
N=1024; 
for t=1:N 
x(t)=bark(ceil(t/Tb)); 
end 
 
% for t=1:20 
%     k=(t-1)*8+1; 
%     if x(k)==1 
%        x(k)=-1;  
%    else 
%        x(k)=1; 
%    end 
% end 
% for t=21:50 
%     k=(t-1)*8+1; 
%     if x(k)==1 
%        x(k)=-1; 
%        x(k+1)=-1; 
%    else 
%        x(k)=1; 
%        x(k+1)=1; 
%    end 
% end 
% for t=51:120 
%     k=(t-1)*8+1; 
%     if x(k)==1 
%        x(k)=-1; 
%        x(k+1)=-1; 
%    else 
%        x(k)=1; 
%        x(k+1)=1; 
%    end 
% end 
%      
snr=10; 
noise=randn(1,N); 
AN=sqrt(std(noise)^2*10^(snr/10)); 
x=AN*x+noise; 
save phase x; 
% [b,a] = butter(5,0.65); 
% x=filter(b,a,x); 
%%%%%%%%%%%%%%%%%%%%%%%%%% 
%求循环谱% 
%%%%%%%%%%%%%%%%%%%%%%%%%% 
alpha_len=(-0.5:1/N:0.5-1/N); 
M=N/16; 
X=fft(x); 
Y=X; 
% figure(5); 
% plot(abs(X)); 
X=fftshift(X); 
figure(4); 
plot(abs(X)); 
 
for alpha=1:N/2-M/2+1, 
    for f1=1:N/2-M/2+2-alpha, 
        if f1==1, 
            tmp(alpha,f1)=X([f1+N/2+alpha-1-M/2:f1+N/2+alpha-2+M/2])*(X([f1+N/2-alpha+1-M/2:f1+N/2-alpha+M/2]))'; 
        else tmp(alpha,f1)=tmp(alpha,f1-1)-X(f1+N/2+alpha-2-M/2)*conj(X(f1+N/2-alpha-M/2))+X(f1+N/2+alpha-2+M/2)*conj(X(f1+N/2-alpha+M/2)); 
        end; 
    end; 
end;%tmp为alpha>0,f>0的四分之一平面 
% for alpha=1:N/2-M/2+1, 
%     for f1=1:N/2-M/2+2-alpha, 
%         if f1==1, 
%             tmp(alpha,f1)=X([f1+N/2+alpha-1-M/2:f1+N/2+alpha-2+M/2])*(X([f1+N/2-alpha+1-M/2:f1+N/2-alpha+M/2]))'; 
%         else tmp(alpha,f1)=tmp(alpha,f1-1)-X(f1+N/2+alpha-2-M/2)*conj(X(f1+N/2-alpha-M/2))+X(f1+N/2+alpha-2+M/2)*conj(X(f1+N/2-alpha+M/2)); 
%         end; 
%     end; 
% end;%tmp为alpha>0,f>0的四分之一平面 
 
S=[tmp([1:N/2-M/2+1],[N/2-M/2+1:-1:2])/M tmp/M];% f轴扩充 
S=[S([N/2-M/2+1:-1:2],:);S]; %alpha轴扩充 
S1=zeros(N,N); 
S1([M/2:N-M/2],[M/2:N-M/2])=S; 
X1=S1(:,N/2+1); 
X2=S1(N/2+1,:); 
S1=abs(S1); 
 
f_axe=[-0.5:1/N:0.5-1/N]*2*pi; 
alpha_axe=(-0.5:1/N:0.5-1/N)*2; 
figure(1); 
% mesh(f_axe,alpha_axe,S1); 
contour(f_axe,alpha_axe,S1); 
xlabel('f'),ylabel('alpha'); 
figure(2) 
plot(alpha_axe,abs(X1)); 
xlabel('alpha') 
figure(3) 
plot(f_axe,abs(X2)); 
xlabel('f') 
 
[temp,J]=max(abs(X1)); 
X1(J)=0; 
[temp,I]=max(abs(X1)); 
xxx=(N/2-I)*2/N 
% J 
% xxx=(J-I)*1/J