clear;close all;clc;

[u,fs,bits]=wavread('D:\matlab\toolbox\vr\vrealm\program\sounds\bird');%讀入語音數據 
N=length(u);

%fs=10000;
t=0:1/fs:(N-1)/fs;%使得t的長度和u的長度一樣，方便時域畫圖
k=0:N-1;
f=fs*k/N;       %使得f的長度和u的長度一樣，方便頻域畫圖
uw=fft(u);     %求出原始信號的頻譜

figure(1);
subplot(2,2,1);plot(t,u,'k');xlabel('t/Sec');ylabel('u(t)');axis([0,0.89,-1,1]);%title('原始信號');;line([0,max(t)],[0,0])
%subplot(2,2,2);plot(f,abs(uw)/max(abs(uw)));xlabel('f/Hz');ylabel('uw(f)');%title('原始信號頻譜');%頻譜歸一化了
ttt=abs(uw)/max(abs(uw));
subplot(2,2,2);plot(f(1:length(f)/2),ttt(1:length(f)/2),'k');xlabel('f/Hz');ylabel('uw(f)');axis([0,5000,0,1]);%title('原始信號頻譜');%頻譜歸一化了
sound(u);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%不同初始條件加密
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xx=zeros(1,N);
yy=zeros(1,N);
xx1=zeros(1,N);
xx21=zeros(1,N);
xx22=zeros(1,N);

xx1(1)=0.655;
xx21(1)=0.2;
xx22(1)=1.2+exp(-11);

for k=2:N
    xx(k)=0.5*(xx1(k-1)-xx(k-1)*xx1(k-1))+0.5*xx21(k-1)*u(k-1);
    yy(k)=xx(k);
    
    xx1(k)=4*xx1(k-1)*xx1(k-1)*xx1(k-1)-3*xx1(k-1);
    xx21(k)=1+0.3*xx22(k-1)-1.4*xx21(k-1)*xx21(k-1);
    xx22(k)=xx21(k-1);
end
yyw=fft(yy);

figure(2);
subplot(2,2,1);plot(t,yy,'k');xlabel('t/Sec');ylabel('yy(t)');axis([0,0.89,-1,1]);%;title('加密信號');axis([0,2,-1,1]);line([0,max(t)],[0,0]);
tt=abs(yyw)/max(abs(yyw));
subplot(2,2,2);plot(f(1:length(f)/2),tt(1:length(f)/2),'k');xlabel('f/Hz');ylabel('yyw(f)');axis([0,5000,0,1]);%title('加密信號頻譜');%頻譜歸一化了
sound(yy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%錯誤解密
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zz1(1)=0.655;
zz21(1)=0.2;
zz22(1)=1.2;

for k=2:N   
    zz1(k)=4*zz1(k-1)*zz1(k-1)*zz1(k-1)-3*zz1(k-1);
    zz21(k)=1+0.3*zz22(k-1)-1.4*zz21(k-1)*zz21(k-1);
    zz22(k)=zz21(k-1);
end

for k=2:N-1
    if abs(zz21(k-1))>0.001
        w2(k)=(2*yy(k+1)-zz1(k)*(1-yy(k)))/zz21(k);
    else 
        w2(k)=w2(k-1);
    end
end
for k=1:N-2
    u2(k)=w2(k+1);%考慮到系統必須是個因果系統，u1(t)即為解密恢復的語音信號
end
u2(N)=0;
uw2=fft(u2);

figure(3);
subplot(2,2,1);plot(t,u2,'k');xlabel('t/Sec');ylabel('u2(t)');axis([0,0.89,-1,1]);%;title('加密信號');axis([0,2,-1,1]);line([0,max(t)],[0,0]);
tt=abs(uw2)/max(abs(uw2));
subplot(2,2,2);plot(f(1:length(f)/2),tt(1:length(f)/2),'k');xlabel('f/Hz');ylabel('uw2(f)');axis([0,5000,0,1]);%title('加密信號頻譜');%頻譜歸一化了
sound(u2);
%wavwrite(u2,'D:\incorrect')