clear

 [ml,fs1,nbits1]=wavread('c:\interim\ml_fft_8k.wav');
 [st,fs2,nbits2]=wavread('c:\interim\st_fir_8k.wav');
 load c:\interim\fir_design_stm32f4\pdmsignal.mat; 
 plot(1:length(st),20*log10(abs(fft(st))),'r.-',1:length(ml),20*log10(abs(fft(ml))),'b.-');
 signal=signal-mean(signal);
 
 %xx=signal(1:127*fix(length(signal)/127));
 %xx=reshape(xx,127,length(xx)/127);
 %xx=sum(xx);
 
 xx=filter(ones(1,512),1,signal);
 
 [fb,fa]=ellip(4,3,100,6000/1e6);
 %fa=1;fb=ones(1,1024);
 xx=filter(fb,fa,[signal signal]);
 %soundsc(xx(1:127:end),fs1);
 soundsc(xx(end/2:128:end),fs1)
 
 H=20*log10(abs(freqz(fb,fa,10000)));
 H=H-max(H);
 plot(H);
 bet=0.9984;wi=1.023; 
 fa1=real(poly([bet*exp(j*(pi/180)*wi) bet*exp(-j*(pi/180)*wi)]));
 bet=1;wi=17.63; 
 fb1=real(poly([bet*exp(j*(pi/180)*wi) bet*exp(-j*(pi/180)*wi)]));
 bet=0.9961;wi=0.426; 
 fa2=real(poly([bet*exp(j*(pi/180)*wi) bet*exp(-j*(pi/180)*wi)]));
 bet=1;wi=7.38; 
 fb2=real(poly([bet*exp(j*(pi/180)*wi) bet*exp(-j*(pi/180)*wi)]));
 %fa1=poly([0.9983+j*0.0178 0.9983-j*0.0178]);
 %fa2=poly([0.9961+j*0.0074 0.9961-j*0.0074]);
 
 fa1(2)=-2+115/32768;
 fa1(3)= 1-105/32768;
 
 fa2(2)=-2+257/32768;
 fa2(3)= 1-255/32768;

 %fb1(2)=-2+3078/32768; 
 fb1(2)=-2+(4096-1024)/32768; 
 fb2(2)=-2+512/32768;
 
 H1=((freqz(fb1,fa1,10000)));
 H2=((freqz(fb2,fa2,10000)));
 HH=H1.*H2;
 HH=20*log10(abs(HH));
 HH=HH-max(HH);
 von=1;bis=3000;
 plot(von:bis,H(von:bis),'b.-',von:bis,HH(von:bis),'r.-')
 
 sisi=[signal signal];
 siso=zeros(size(sisi));
 for(k=3:length(sisi))
   siso(k)=        fb1(1)*sisi(k);
   siso(k)=siso(k)+fb1(2)*sisi(k-1);
   siso(k)=siso(k)+fb1(3)*sisi(k-2);
   siso(k)=siso(k)-fa1(2)*siso(k-1);
   siso(k)=siso(k)-fa1(3)*siso(k-2);  
 end
 sisi=siso;
 for(k=3:length(sisi))
   siso(k)=        fb2(1)*sisi(k);
   siso(k)=siso(k)+fb2(2)*sisi(k-1);
   siso(k)=siso(k)+fb2(3)*sisi(k-2);
   siso(k)=siso(k)-fa2(2)*siso(k-1);
   siso(k)=siso(k)-fa2(3)*siso(k-2);  
 end
 
 siso=siso(100000:end);
 siso=siso/max(abs(siso));
 
 soundsc(siso(1:128:end),fs1)
 
 return
 