afrplot.m

binDCT fron D.Tran used in Image compression

function [DC, smax, senergy] = afrplot(P,G, odd, plotornot)
% Function to plot the frequency response of the analysis FB as well as  
% its coding gain, DC attenuation, frequency attenuation, and mirror
% frequency attenuation

% Trac D. Tran 
% ECE Department, The Johns Hopkins University
% 3400 North Charles Street, 105 Barton Hall,
% Baltimore, MD 21218
% E-mail: [e-mail]
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Copyright (c) 2000 Trac D Tran
% This program is Copyright (c) by Trac D Tran.
% It may not be redistributed without the consent of the copyright
% holders. In no circumstances may the copyright notice be removed.
% The program may not be sold for profit nor may they be incorporated
% in commercial programs without the written permission of the copyright
% holders. This program is provided as is, without any express or
% implied warranty, without even the warranty of fitness for a
% particular purpose.
%-----------------------------------------------------------------------

if nargin == 2
    odd = 0;
    plotornot = 1;
end    

if nargin == 3
    plotornot = 1;
end    

[M,L]=size(P);
F=abs(fft(P',512));
F=(F+(F<1e-20)*1e-20)';
Forg = F(:, 1:256) / F(1,1); % no log
F=20*log10(F(:,1:256)/F(1,1));
DC=max(F(2:M,1));
%F=max(-50,F);
N=256/M;bis=2*N+1;

[tmp,imax]=max(F');
Fdiff=F-[F(:,1) F(:,1:255)];smax=-50;
for i=1:M
   iright=find(Fdiff(i,imax(i):256)>3);
   if ~isempty(iright) %iright~=[]
      iright=iright(1);
      smax=max(smax,max(F(i,imax(i)+iright:256)));
   end;
   ileft=find(fliplr(Fdiff(i,1:imax(i)))<-3);
   if ~isempty(ileft) %ileft~=[]
      ileft=ileft(1);
      smax=max(smax,max(F(i,1:imax(i)-ileft)));
   end;
end;

if plotornot == 1
    x=0:1/512:.5-1/512;
    figure;
    plot(x,F');
    axis([0 0.5 -40 5]);
%    axis([0 0.5 -100 5]);
    set(gca,'fontsize',7);

    if odd == 0
        title(['DC Att. \geq ',num2str(-DC),', Mirr Att. \geq ',...
        num2str(-mirrfreq(P)),', Stopband \geq ',num2str(-smax),...
        ', Coding Gain = ',num2str(bgtc(P,G)),' dB']);
    else
        %don't use mirror freq. for odd channel
        title(['DC Att. \geq ',num2str(-DC), ', Stopband \geq ',num2str(-smax),...
        ', Coding Gain = ',num2str(bgtc(P,G)),' dB']);
    end

    tt=get(gca,'title');set(tt,'fontsize',10);
    xlabel('Normalized Frequency');
    xl=get(gca,'xlabel');set(xl,'fontsize',11);
    ylabel('Magnitude Response (dB)');
    yl=get(gca,'ylabel');set(yl,'fontsize',11);
    grid on;
end

senergy = smax;

break;

% stopband energy
%find -3db point.
e_stopband = zeros(1, M);
hfbandwth = round(1/(2*(M-1)) * 256);
for idx = 1 : M
    sbctr = round((idx - 1) / (M-1) * 256);
    if idx == 1
        sum_tail = 256;
        stopbandctr = sbctr + hfbandwth;
        window_head = max(1, stopbandctr - hfbandwth);
        window_tail = min(256, stopbandctr + hfbandwth);
        tmp = abs(Forg(idx, window_head : window_tail)  - 0.7071);     % min value is at the -3db point.
        sum_head = (window_head - 1) + find(tmp == min(tmp));
        e_stopband(idx) = sum( 10*log10(Forg(idx, sum_head : sum_tail) .^ 2)) * pi / 256;
    elseif idx == M
        sum_head = 1;
        stopbandctr = sbctr - hfbandwth;
        window_head = max(1, stopbandctr - hfbandwth);
        window_tail = min(256, stopbandctr + hfbandwth);
        tmp = abs(Forg(idx, window_head : window_tail)  - 0.7071);
        sum_tail = (window_head - 1) + find(tmp == min(tmp));
        e_stopband(idx) = sum( 10*log10(Forg(idx, sum_head : sum_tail) .^ 2)) * pi / 256;
    else
        %low freq side
        sum_head = 1;
        stopbandctr = sbctr - hfbandwth;
        window_head = max(1, stopbandctr - hfbandwth);
        window_tail = min(256, stopbandctr + hfbandwth);
        tmp = abs(Forg(idx, window_head : window_tail)  - 0.7071);
        sum_tail = (window_head - 1) + find(tmp == min(tmp));
        e_stopband(idx) = sum( 10*log10(Forg(idx, sum_head : sum_tail) .^ 2)) * pi / 256;

        %high freq side
        sum_tail = 256;
        stopbandctr = sbctr + hfbandwth;
        window_head = max(1, stopbandctr - hfbandwth);
        window_tail = min(256, stopbandctr + hfbandwth);
        tmp = abs(Forg(idx, window_head : window_tail)  - 0.7071);     % min value is at the -3db point.
        sum_head = (window_head - 1) + find(tmp == min(tmp));
        e_stopband(idx) = e_stopband(idx) + sum( 10*log10(Forg(idx, sum_head : sum_tail) .^ 2)) * pi / 256;
    end
end

senergy = sum(e_stopband);
if plotornot == 1
    s = sprintf('Stopband energy: %f.', senergy);
    disp(s);
end