% Verifica della energy-compaction property della DCT-2
% usando lo stesso segnale del libro Oppenheim-Schafer

clear all;

% generazione del segnale
N= 32;	% durata segnale: DEVE ESSERE PARI, PER IL PROGRAMMA (v. dopo)
a= 0.9; w0= 0.1*pi;	% parametri del segnale
nvec= (0:N-1);	% ascisse
x= (a.^nvec).*cos(w0*nvec);	% sinusoide a decadimento esponenziale

% visualizzazione del segnale
figure(1); clf;
stem(nvec,x,'m');
xlabel('n'); ylabel('x(n)');
title('segnale');
axis([0 N-1 -1 1]);

% visualizzazione dello spettro DFT e DCT-2
Xk= fft(x);
Xc2= dct(x);
% attenzione: ? definita come DCT-2 Unitaria
% i valori della unitaria sono pari a quelli che conosciamo, ma
% divisi per sqrt(2*N), tranne il primo, che ? diviso per 2*sqrt(N)
figure(2); clf;
subplot(3,1,1);
stem(nvec,real(Xk),'b');
xlabel('k'); ylabel('Re[X(k)]  (DFT)');
axis([0 N/2 0 6]);
subplot(3,1,2);
stem(nvec,imag(Xk),'b');
xlabel('k'); ylabel('Im[X(k)]  (DFT)');
axis([0 N/2 -10 5]);
subplot(3,1,3);
stem(nvec,Xc2,'r');
xlabel('k'); ylabel('X^{c2}(k)  (DCT)');
axis([0 N-1 -10/sqrt(2*N) 20/sqrt(2*N)]);

% calcolo l'errore Mean Square di approssimazione DFT, al variare del numero m di coeff. azzerati
ErrDFT= zeros(1,N/2); % preparazione vettore valori errore
for m=1:2:N-1,
  XkT= [Xk(1:N/2-(m+1)/2) zeros(1,m) Xk(N/2+(m+1)/2:N)];
  xT= real(ifft(XkT));	% real per evitare errori di approssimazione
  ErrDFT((m+1)/2)= sum(abs(x-xT).^2)/N;	% errore mean square
end;	% for m
mvec= 1:2:N-1;	% vettore ascisse per ErrDFT

% calcolo l'errore Mean Square di approssimazione DCT-2, al variare del numero m di coeff. azzerati
ErrDCT= zeros(1,N-1); % preparazione vettore valori errore
for m=1:N-1,
  Xc2T= [Xc2(1:N-m) zeros(1,m)];
  xTc2= real(idct(Xc2T));	% real per evitare errori di approssimazione
  ErrDCT(m)= sum(abs(x-xTc2).^2)/N;	% errore mean square
end;	% for m
mvecc2= 1:N-1;	% vettore ascisse per ErrDCT


% confronto errori
figure(3); clf;
plot(mvecc2,ErrDCT,'r.-');
hold on;
plot(mvec,ErrDFT,'bo--');
xlabel('n. coefficienti posti a zero (m)'); ylabel('Errore MS');
title('Energy Compaction della DCT-2');
legend('DCT-2 errore di troncamento','DFT errore di troncamento',2);
axis([0 N 0 0.1]);