
clc;clear all;close all;
[x,Fs,nbits]=WAVREAD('test1.wav');
[xx,Fs,nbits]=WAVREAD('Bill2.wav');
wstring='WATERMARK!';
temp=strcat(wstring,'-----------');        
watermark(1:10)=temp(1:10);

qq=dec2bin(watermark,8);
bits=[1 1 1 1 1 1 1 1 ];
for K=1:size(qq,1),
    for L=1:8,
        bits(L+(K-1)*8+8)=str2num(qq(K,L)); %
    end
end
  
delay0=128;
delay1=256;
nbits=88;
N=512;
T=0.5;
Nram = 1+floor((size(x,1)-N/2)/N); %pocet ramcu
q=ceil(size(x,1)/N)-1;
matrix=zeros(N,Nram);
matrixWTM=zeros(N,Nram);
energy=zeros(q,1);
zer=zeros(nbits,1);
one=zeros(nbits,1);
energy2=zeros(q,1);
y=zeros(N+(Nram*N),2);
x0=zeros(1,Nram);
x1=zeros(1,Nram);
ffty = zeros(N,1);
ifftlogy = zeros(1,N);
ifftlogy1 = zeros(1,N);
sync=zeros(1,8*N);
recbits = zeros(1,Nram);
xcorrel0 = zeros(1,N);
xcorrel1 = zeros(1,2*N);
n=1;
counter =1;
startblok=0;
stopblok=0;
figure(1);


fs1=512;
fs2=1/512;
Q0=0.5;
Q1=0.5;
fb1=fs1/Q0;
fb2=fs2/Q1;
Fvz=44100;
NN=2048;
c=(tan(pi*fb1/Fvz)-1)/(tan(pi*fb1/Fvz)+1);
d=-cos(2*pi*fs1/Fvz);
Ab=[-c,d*(1-c),1];
Aa=[1,d*(1-c),-c];
[Hb1,Ha1]=parallel([1 0],[1 0],-Ab,Aa);
Hb1=0.5*Hb1;
Ha1=0.5*Ha1;




%--- tvorba mikrosegmentų ---
%matrix = ramce(x,512,256);
%[q,r]= size(matrix);
%okno = hamming(q)*ones(1,r);
%matrix = matrix.*okno;

L1=10;
L2=502;
L3=N;



       
%synchronizace
synch=zeros(N,1);
for j=1:N-1
 for i = 1:q
     if (x(j+N*(i-1)) < 0 & x(j+1+N*(i-1)) > 0)
         synch(j)= synch(j)+1;
     elseif (x(j+N*(i-1)) > 0 & x(j+1+N*(i-1)) < 0)
         synch(j)= synch(j)+1;
     end
end
end
syncindex=0;
s_tmp=0;
for j=1:N
    if(synch(j)> s_tmp)
    s_tmp=synch(j);
    syncindex=j;
    end
end
    
   
d=0;%syncindex;
for i=1:(q-1)
 for j = 1:N
  matrix(j,i) = x(j+d+N*(i-1));
  %matrix2(j,i) = x2(j+d+N*(i-1));
  energy(i)=energy(i)+ matrix(j,i)^2;
  %energy2(i)=energy2(i)+ matrix2(j,i)^2;
 end
   energy(i)=energy(i)/N;
   %energy2(i)=energy2(i)/N;
   if(energy(i)>1e-5) & (startblok==0)
       startblok=i+7;
   end
   if(energy(i)<1e-5) & (startblok~=0) | (stopblok==0) & (i>q-10)
   stopblok=i-1;
   end
end
a=1;
t=linspace(0,N*1000/Fvz,N);


for i=startblok:stopblok %sloupce

        ffty=log(abs(fft(matrix(:,i),N).^2));
        ifftlogy = real(ifft(ffty));
        %fftyWTM=log(abs(fft(matrixWTM(:,i),N).^2));
        %ifftlogyWTM = real(ifft(fftyWTM));

        xcorrel0 = abs(xcorr2(ifftlogy));
        xcorrel0(N-2:N+2)=0;

        
        %figure(1);
       % subplot(311);
       % plot(t,xcorrel0(1:N));axis xy;axis tight;xlabel('Time[ms]');ylabel('Amplitude');
        %line([delay0 delay0],[0.5 5]);
        %line([delay1 delay1],[0.5 5]);
 
        
        a=a+1;
        %pause;
        
        str=0;
        ee=0;
    x0(counter)=sum(xcorrel0(str+delay0-ee:1:str+delay0+ee));
    x1(counter)=sum(xcorrel0(str+delay1-ee:1:str+delay1+ee));
   
    if(x0(counter)<x1(counter))
        
        recbits(1,counter)=1;
   
    end
    counter = counter +1;
end


for i = 1:floor(stopblok/nbits)
   for j=1:nbits
if(recbits((i-1)*nbits+j)==1)
    if(x1((i-1)*nbits+j)>T)
    one(j) = one(j)+1;
    end
    
else    
    if(x0((i-1)*nbits+j)>T)
        zer(j) = zer(j)+1;
    end
end
   end
end
xxxx=sum(one(:,1))+sum(zer(:,1));
watermark = zeros(size(zer,1),1);
for j=1:size(zer,1)
if(one(j)>zer(j))
watermark(j)=1;
end
end
recbits
watermark'

%SNR

SNR = 10*log10(((sum(x(:,1).^2)-sum(xx(:,1).^2))./sum(xx(:,1).^2)))
%BER%%%%%%%%%%%
BER=0;
for i = 1:floor(counter/nbits)
   for j=1:nbits
    if(recbits((i-1)*nbits+j)~=bits(j))
    BER=BER+1;
    end
   end
end
BER=(BER/(floor(q/nbits)*nbits))*100
