function psycho
%Test_MPEG
%
%-------------------------------------------------------------------------------
common;

% Build a demo sample
fs = 44100;
f = [1500 1700 1800 2000 10000 11000 12000 13000];
a = [   1    2    1   .5     2   1.5     3     1];
t = 0:1/fs:1;
x = randn(size(t));
for i = 1:length(f),
   x = x + a(i) * sin(2 * pi * f(i) * t);
end
x = x / max(abs(x));
%[x fs] = wavread('svega');

% Load tables.
[TH, Map, LTq] = absolth(1, fs, 128); % Threshold in quiet
CB = critband(1, fs);
C = anawind;

% Process the input vector x.
for OFFSET = 1:384:length(x),
   S = [];

%%% Subband filter analysis. Layer 1 uses 12 samples per subband.
   
   % Analysis subband filtering [1, pp. 67].
   for i = 0:11,
      S = [S; ansbbfil(x, OFFSET + 32 * i, C)];
   end
   
   % Scalefactor calculation [1, pp. 70].
   scf = scalefc(S);
   
%%% Psychoacoustic analysis.

	% Compute the FFT for time frequency conversion [1, pp. 110].
	X = fftana(x, OFFSET);
   
   % Determine the sound pressure level in each  subband [1, pp. 110].
   Lsb = soundpr(X, scf);
   
   % Find the tonal (sine like) and non-tonal (noise like) components
   % of the signal [1, pp. 111--113]
   [Flags Tonal_list Non_tonal_list] = tonal(X, TH, Map, CB);
   
   % Decimate the maskers: eliminate all irrelevant maskers [1, pp. 114]
   [Flags Tonal_list Non_tonal_list] = ...
      decimat(X, Tonal_list, Non_tonal_list, Flags, TH, Map);
   
   % Compute the individual masking thresholds [1, pp. 113--114]
   [LTt, LTn] = ...
      indivith(X, Tonal_list, Non_tonal_list, TH, Map);
   
   % Compute the global masking threshold [1, pp. 114]
   LTg = globalth(LTq, LTt, LTn);
   
   if (DRAW)
      disp('Global masking threshold');
      hold on;
      plot(TH(:, INDEX), LTg, 'k--');
      hold off;
      title('Masking components and masking thresholds.');
   end
   
   % Determine the minimum masking threshold in each subband [1, pp. 114]
   LTmin = minth(LTg, Map);
   if (DRAW)
      figure; plot(LTmin); title('Minimum masking threshold');
      xlabel('Subband number'); ylabel('dB'); pause;
   end
   
   % Compute the singal-to-maks ratio
   SMRsb = Lsb - LTmin;
end