MP3原理和实现代码3

function [SMR min_threshold_subband frame_psd_dBSPL ...
    masking_threshold] = MPEG1_psycho_acoustic_model1（frame）

% function [SMR min_threshold_subband frame_psd_SPL ...
%    masking_threshold] = MPEG1_psycho_acoustic_model1（frame）
% Computes the masking threshold （in dB） corresponding to psycho-acoustic
% model #1 used in MPEG-1 Audio （cf  ISO/CEI  norm 11172-3:1993 （F） pp.
% 122-128）. 
% Input |frame| length should be 512 samples in the [-1+1] range. 
% |SMR| returns 27 signal-to-mask ratios （in dB）.
% |min_threshold_subband| returns the minimun of |masking threshold| in each
% of the 32 sub-bands. 
% |frame_psd_SPL| returns the estimated PSD of the input frame in dB SPL
% assuming the level of full scale signals is set to 96 dB SPL.
%
% Copyright N. Moreau ENST Paris 19/03/02
% Modified by Thierry Dutoit FPMs Mons 03/05/07

global LTq_i LTq_k Table_z Frontieres_i Frontieres_k Larg_f

if size（LTq_i）==[00]
   MPEG1_psycho_acoustic_model1_init;
end;

N = length（frame）;
if N ~= 512
    disp（‘frame length must be set to 512‘）
    return
end

% FFT
% ***

hann = sqrt（8/3）/2*[ones（N 1） - cos（2*pi*（0:N-1）‘/N）];
if sum（abs（frame）） > 0
    X1 = fft（frame.*hann）;
    X1 = （abs（X1（1:N/2+1））.^2）/N;
    perio_xn_db = 10*log10（X1）;
else
    perio_xn_db = zeros（N/2+11）;
end
%offset = max（perio_xn_db） - 96;
%X = perio_xn_db - offset;

% NB: since the absolute acoustic level set by the listener is not known by
% the MPEG psycho-acoustic model it assumes that the level is set such
% that  a full-scale signal corresponds to 96 dB SPL. Since 16 bits signals
% have about 96 dB of dynamics this implies that the LSB is close to 
% the absolute auditory threshold.
% Since the absolute value of input samples is assumed to be <1 a
% full-scale signal i.e. ones（1:512） will produce a PSD peak at
% 10*log10（512）=27.09dB. Hence the 96-27.09 dB offset.
offset=96-27.09;
X = perio_xn_db + offset;
frame_psd_dBSPL=X（1:256）;

% Tonal and noise masker detection 
% *************************

% Local maximum search 

max_local = zeros（250 1）;
for k = 3:250
    if X（k） > X（k-1） & X（k） >= X（k+1）
        max_local（k） = 1;
    end
end

tonal = zeros（250 1）;
for k = 3:62
    if max_local（k）
        tonal（k） = 1;
        for j = [-2 2]
            if X（k） - X（k+j） < 7
                tonal（k） = 0;
            end
        end
    end
end

for k = 63:126
    if max_local（k）
        tonal（k） = 1;
        for j = [-3 -2 2 3]
            if X（k） - X（k+j） < 7
                tonal（k） = 0;
            end
        end
    end
end

for k = 127:250
    if max_local（k）
        tonal（k） = 1;
        for j = [-6:-2 2:6]
            if X（k） - X（k+j） < 7
                tonal（k） = 0;
            end
        end
    end	
end

% Tonal masker detection

X_tm = zeros（2501）;

for k = 1:250
    if tonal（k）
        temp = 10^（X（k-1）/10） + 10^（X（k）/10） + 10^（X（k+1）/10）;