cyclic2.m

基于matlab实现的循环谱估计算法

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