CZS cvičení 9(Matlab)

% cvičení 9

% Nacteni/generovani signalu

% Generovani sinusovky
f=100;
fs=8000;
tmax=0.4 ;
%sig=sin(2*pi*f*(0:1/fs:tmax)) ;

% Nacteni cisteho recoveho signalu
% sig=loadbin('vm0.bin');  fs=16000 ;
% sig=loadbin('SA110992.CS0');  fs=16000 ;
% sig=loadbin('ma034014.ils');  fs=16000 ;



% Nacteni recoveho signalu se sumem
 sig=loadbin('SA110992_auto1.CS0');  fs=16000 ;
% sig=loadbin('SA107S06_auto2.CS0');  fs=16000 ;
% sig=loadbin('ma034014_auto3.ils');  fs=16000 ;



% Vypocet poctu oken
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sig=sig(:);
slen=length(sig);

wlen = 256;
wstep = 128;
wnum = (slen-wlen)/wstep+1;


w=hamming(wlen);
% w=rectwin(wlen);
% w=hanning(wlen);
% w=blackman(wlen);

% Vykresleni vstupniho signalu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); clf ;
subplot(211);
plot(sig);
title('Vstupni signal');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vytvoreni masky pamove propusti
%   fc1 = 300 Hz
%   fc2 = 3400 Hz

NDFT=512;

fc1 = 300 ;
fc2 = 300 ;%3400 ;

k1=floor((fc1/fs)*NDFT);
k2=floor((fc2/fs)*NDFT);

k=k2-k1;
H=[[zeros(1,k1)] [ones(1,k2-k1)] [zeros(1,(NDFT-2*k2)+1)] [ones(1,k2-k1)] [zeros(1,k1-1)]];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nalezeni odhadu amplitudoveho spektra sumu
%   z prvniho casti signalu bez reci
%   (prvnich 30 segmentu)





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generovani nulove vystupni posloupnosti
out = zeros(slen,1);

winit = 30;

Nabs = [];

% toto je welchova metoda
for i=1:winit;
ii=(i-1)*wstep+1;
  jj=(i-1)*wstep+wlen;

  frame=sig(ii:jj).*w;
 
  X=fft(frame);
  Xabs = abs(X) ;
  Xangle = angle(X) ;
  Nabs = [Nabs Xabs];
end

Npsd = pwelch(sig(1:500),w,[],[],'twosided');
%Gx = abs()/n;

Navg = mean(Nabs,2);
 


% Hlavni cyklus
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:wnum,

  % Nastaveni indexu a vyber segmentu
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ii=(i-1)*wstep+1;
  jj=(i-1)*wstep+wlen;

  frame=sig(ii:jj).*w;
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Zde bude modifikace
  X=fft(frame);
  Xabs = abs(X) ;
  Xangle = angle(X) ;
  Nabs = [Nabs Xabs];
 
  % 1) Odladeni bez modifikace

  %yframe = frame ;
 
  Yabs = Xabs - Navg;
  Yabs = (Yabs+ abs(Yabs))./2;
  Yabs = abs(Yabs);
 
  yframe = real(ifft(Yabs.*exp(j*Xangle)));
  % 2) pasmova filtrace
 
 
  % 3) spektralni odecitani
 

 
  % Pricteni do vystupni posloupnosti (secteni presahu)
  out(ii:jj)=out(ii:jj)+yframe ;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Normovani vystupniho signalu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% scale = ....
% out = out ./ ..... ;


subplot(212);
plot(out) ;
title('Vystupni signal');

figure(2)
subplot(211);
spectrogram(sig,wlen,[],[],fs,'yaxis');
title('Spektrogram vstupniho signalu');
subplot(212);
spectrogram(out,wlen,[],[],fs,'yaxis');
title('Spektrogram vystupniho signalu');
%sound(sig,fs);
sound(out,fs);