clc,clear
[x, fs]=wavread('no1.wav');
sound(x, fs);			              % 播放此音頻
time=(1:1024)/fs;	                  % 時間軸的向量
figure(1);
plotspec(x,1/fs);			          % 時域和頻域波形

%============================================================

N=128;                                %濾波器的長度
NUM=1000;
Ts=1/fs;
%Ts=1/NUM;
div=0.1;                              %截止頻率的初始步長
eval=exp(-j*pi*(0:N-1)/2);            %在ω=n*pi/2處的復指數值        
freq=0.5;                             %初始化截止頻率
target=1/sqrt(2);                     %在ω=n*pi/2處,希望得到的幅值

h0=fir1(N-1,freq,blackman(N));        %用窗函數法找到濾波器系數
figure(2);plotspec(h0,Ts);

 lenH=length(h0);
 h1=zeros(1,lenH);
 for n=1:lenH
     h1(n)=(-1)^(n-1)*h0(n);
 end
 
   figure(3);plotspec(h1,Ts);         % 濾波器的時域和頻域波形
   x0=filter(h0,1,x);
   x1=filter(h1,1,x);
   figure(4); plotspec(x0,Ts);        %第一路濾波信號的時域和頻域波形
   figure(5); plotspec(x1,Ts);        %第二路濾波信號的時域和頻域波形
   
g0=zeros(1,lenH);
g1=zeros(1,lenH);
for n=1:lenH
     g0(n)=2*h0(n);
     g1(n)=-2*((-1)^(n-1))*h0(n);
end

[x00,r1] = DPCM_function(x0);              %DPCM編碼
[x11,r2] = DPCM_function(x1);              %DPCM編碼


figure(6); plotspec(x00,Ts);       %第一路信號DPCM編碼后的時域和頻域波形
figure(7); plotspec(x11,Ts);       %第二路信號DPCM編碼后的時域和頻域波形
y0=filter(g0,1,x00);
y1=filter(g1,1,x11);
figure(8); plotspec(y0,Ts);           %第一路濾波信號的時域和頻域波形
figure(9); plotspec(y1,Ts);           %第二路濾波信號的時域和頻域波形
X=y0+y1;
figure(10); plotspec(X,Ts);           %最后處理后的時域和頻域波形
sound(X, fs);			              %播放此音頻

%=====================求編碼速率=============================================
r1
r2