% clear all
%水印嵌入程序
%通過奇偶嵌入水印
%讀取掩護媒體
audio_name='D:\sound\鋼琴曲 - 少女的祈禱.wav';
[y,fs,br]=wavread(audio_name);
%若掩護媒體為雙聲道則取第一個聲道,在audio中嵌入水印
audio=y(:,1);
audio_length=size(audio,1);%掩護媒體長度
audio_new=audio;%audio_new為嵌入水印的音頻媒體


image_read;%讀取水印圖象


sect=1024;%將掩護媒體分幀的每幀長度為1024
k=100000;%起始點
w=1;%w為水印標號
N=1;%N為每幀的嵌入水印個數
s=0.025;%s為嵌入強度系數
%h為漢明窗
i=1;
for i=0:1023
    h(i+1)=sqrt(8/3)*0.5*(1-cos(2*pi*i/1024));
end

while k<=audio_length&w<=water_mark_length
    for i=1:sect
        audio1(i)=audio(k,1);
        k=k+1;
    end
    power=sum(audio1.^2);
    %只在能量>0。02的幀中嵌入水印
    if power>0.02
    %對每幀進行分析，判斷水印嵌入位數
    audio2=audio1.*h;
    f=abs(fft(audio2)).^2;
    t=10*log10(geomean(f(1:128))/mean(f(1:128)));
    if t>=-14.5
        N=4;
    else
        N=0;
        k=k-512;
    end
    if N>0
    for i=1:N
        if w<=water_mark_length
        water_mark1(i)=water_mark(w);
        w=w+1;
        end
    end
%     N=length(water_mark1);

%三級小波分解
    [C,L]=wavedec(audio1,3,'db1');
    a1=appcoef(C,L,'db1',3);%a1為小波三級分解的細節分量
 
    a2=qianru_3(a1,N,water_mark1,s);
    for i=1:L(1)
        C(i)=a2(i);
    end
    iaudio1=waverec(C,L,'db1');
    xx=1;

    for iii=k-sect:k-1
        audio_new(iii)=iaudio1(xx);
        xx=xx+1;
    end
    end
    end

    
    
end
%嵌入水印的音頻文件名
audio_new_name='D:\sound\鋼琴曲 - 少女的祈禱_new22.wav';
wavwrite(audio_new(1:k+40000),fs,audio_new_name);

%畫時域圖
% figure(1)
% subplot(2,1,1)
% plot(audio(70000:k))
% title('原始信號')
% subplot(2,1,2)
% plot(audio_new(70000:k))
% title('嵌入水印信號')
% subplot(3,1,3)
% plot(audio(70000:k)-audio_new(70000:k))
% title('誤差')
%畫頻譜圖
% cc=1:k/2;
% fr=(0:k/2-1)*2/k*fs;
% f1=fft(audio(1:k));
% figure(1)
% subplot(2,1,1)
% plot(fr,abs(f1(cc)))
% title('原始信號頻譜')
% % figure(2)
% subplot(2,1,2)
% f2=fft(audio_new(1:k));
% plot(fr,abs(f2(cc)))
% title('嵌入水印信號的頻譜')
% figure(3)
% plot(fr,abs(f2(cc)-f1(cc)))
% ylabel('誤差');
% xlabel('頻率/HZ');
% title('頻譜誤差')