% $Id: EQ.m 1.9 2018/10/30 03:44:33 rc902779 Exp rc902779 $

% This Matlab program lets the user do the following:

% (1) control the gains, center frequencies, and Q factors of 10 biquad

%     digital filters in a 10-band parametric equalizer (EQ),

% (2) plot the composite frequency response of the 10-band EQ, and

% (3) write the 10 biquad filter coefficients to a header file with the C

%     language syntax so a C program can include the *.h file to implement

%     a fixed 10-band parametric EQ designed by this Matlab program.

%

% Copyright (c) 2018 by Broadcom Inc., All Rights Reserved.



SF = 8000;                  % Sampling Frequency in Hz

%SF = 16000;                 % Sampling Frequency in Hz

N_Stages = 10;              % Number of biquad filter STAGES

Default_Q = 1.5;            % DEFAULT initial Q factor for all biquad filters

Gain = zeros(N_Stages,1);   % initialize array of biquad filter GAINs to 0 dB

%Gain = [-10,3,4,3,-10,-10,-9,11,6,-5]'; % a random example of biquad filter GAINs

%Gain = [-2,1,1,1.5,-1,-2,-1,3,2,-2]'; % a random example of biquad filter GAINs

%Gain = [2,-1,-1,-1.5,1,2,1,-3,-2,2]'; % a random example of biquad filter GAINs

%Gain = [20,0,-2,-15,-20,-10,10,5,-3,-6]'; % a random example of biquad filter GAINs

%Gain = [5,1,-3,-3,-6,-3,7,6,5,-15]'; % a random example of biquad filter GAINs

Q = Default_Q*ones(N_Stages,1);  % array of initial defualt Q factors

PRESELECT_CENTER_FREQ = 1;  % = 1 if assign pre-selected center frequencies

                            % = 0 to let program calculate center frequencies



% Pre-select or calculate the Center_Freq() array

if PRESELECT_CENTER_FREQ == 1,

   if SF == 16000,

      Center_Freq = [300 500 800 1000 2000 3000 4000 5000 6000 7000]';

   elseif SF == 8000,

      Center_Freq = [300 500 800 1000 1500 2000 2500 3000 3500 4000]';

   end;

else

   % Now calculate the default initial center frequencies by dividing the

   % logarithmic frequency range from 20 Hz to 20 kHz (or half the sampling

   % frequency, whichever is smaller) into N_Stages equal frequency bands and

   % then getting the center frequency of each band.

   Min_Freq = 20;              % MINimum FREQuency (Hz) of human hearing range

   Max_Freq = min(20000,SF/2); % MAXimum FREQuency (Hz) of hearing or SF/2

   Min_Log_Freq = log2(Min_Freq); % MINimum LOGarithmic FREQuency

   Max_Log_Freq = log2(Max_Freq); % MINimum LOGarithmic FREQuency

   Log_Freq_BW = (Max_Log_Freq - Min_Log_Freq)/N_Stages; % LOG FREQ BandWidth

   Log_Center_Freq=((Min_Log_Freq+(Log_Freq_BW/2)):Log_Freq_BW:Max_Log_Freq)';

                                          % LOGarithmic CENTER FREQuencies

   Center_Freq = 2.^Log_Center_Freq; % CENTER FREQuencies in Hz

end;



% Convert center frequencies from Hz to normalized frequency (1 = Nyquist freq)

Nyquist_Freq = SF/2;    % NYQUIST FREQuency = half the sampling frequency

N_Center_Freq = Center_Freq/Nyquist_Freq;



% Convert Q factors to bandwidths

Bandwidth = N_Center_Freq ./ Q;



% Declare the biquad filter coefficient arrays

Biquad_Coeff_a = zeros(N_Stages,3);  % for all-pole portion of the filter

Biquad_Coeff_b = zeros(N_Stages,3);  % for all-zero portion of the filter



% Do some preparation for plotting

N_Freq_Points = 16384; % Number of FREQuency POINTS to plot frequency responses

h=zeros(N_Freq_Points,N_Stages); % initialize response of individual biquads

h_cascaded=zeros(N_Freq_Points,1); % initialize response of cascaded filters

V_Max = 0;             % initialize Vertical plotting range MAXimum

V_Min = 0;             % initialize Vertical plotting range MINimum



% Prepare for fixed-point conversion of filter coefficients and gains

Gain_Shift = zeros(N_Stages,1); % # of bits for filter GAIN to SHIFT to get Q15

Q30Scale = 2^30;    % Scale factor to scale floating-point to Q30 format

Q15Scale = 2^15;    % Scale factor to scale floating-point to Q15 format

gain = int16(0);    % Gain shift amount applied after the last section



% Loop through each band and design the biquad filter for that band

for k=1:N_Stages,



   % Call Matlab function for IIR PARAMetric EQ design

   [SOS,SV] = iirparameq(2,Gain(k),N_Center_Freq(k),Bandwidth(k));



   % Convert coefficients of designed biquad filter to 32-bit fixed point Q30

   coef_sos(k,:) = int32(round(Q30Scale*[SOS(1:3) SOS(5:6)])); % Q30 filter coefficients



   % Convert the filter gain to Q15 and record the number of bits of shift

   N_Integer_Bits = 1 + ceil(log2(SV(1))); % Number of INTEGER BITS (w/ sign)

   Q_Format = 16 - N_Integer_Bits; % Q FORMAT (number of fractional bits)

   coef_g(k) = int16(round(((2^Q_Format)-1)*SV(1))); % filter gain fully normalized

   Gain_Shift(k) = 15 - Q_Format;  % record # of bits gain shifted right for

                                   % coef_g(k) to be regarded as a Q15 number

   gain = gain + Gain_Shift(k);    % update gain shift amount after last stage



   % Convert fixed-point filter parameters back to floating-point for plotting

   Biquad_Coeff_a(k,1:3) = [1 double(coef_sos(k,4:5))/Q30Scale];% assign all-pole portion

   g = (2^Gain_Shift(k))*(double(coef_g(k))/Q15Scale); % gain of this stage

   Biquad_Coeff_b(k,1:3) = double(g)*double(coef_sos(k,1:3))/Q30Scale; % assign all-zero portion



   % Compute frequency response of this biquad filter

   [H,w]=freqz(Biquad_Coeff_b(k,1:3),Biquad_Coeff_a(k,1:3),N_Freq_Points);

   h(:,k)=20*log10(abs(H));  % calculate biquad magnitude response in dB

   h_cascaded = h_cascaded + h(:,k); % response of cascaded biquads so far



   % Update the Maximum and Minimum of the frequency responses

   V_Max = max([V_Max;max([h(:,k); h_cascaded])]); % update Vertical MAXimum

   V_Min = min([V_Min;min([h(:,k); h_cascaded])]); % update Vertical MINimum



end;



% Check to make sure that the gain shift amount is within allowed range

if gain > 15,

    fprintf('Gain shift amount exceeds maximum allowed 15 bits!\n');

    fprintf('Please reduce the gain of the biquad filters and try again.\n');

    return;

end;

if gain < -48,

    fprintf('Gain shift amount falls below minimum allowed -48 bits!\n');

    fprintf('Please increase the gain of the biquad filters and try again.\n');

    return;

end;



% Now plot the frequency responses of individual biquad filters (in green)

% and the frequency response of the overall cascaded filter (in red)

figure(1);           % choose Figure 1 for plotting

clf;                 % clear the figure

H_Min = 20;           % set Horizontal plotting range MINimum to 1 Hz

V_Max = V_Max + 1;   % add 1 dB of vertical headroom in plotting

V_Min = V_Min - 1;   % add 1 dB of vertical headroom in plotting

f=w*Nyquist_Freq/pi; % convert angular frequency to linear frequency

semilogx(f,h,'g--',f,h_cascaded,'r'); % plot frequency responses

axis([H_Min Nyquist_Freq V_Min V_Max]); % set the plotting range

xlabel('Frequency (Hz)');

ylabel('Magnitude (dB)');

title('Frequency Responses of Cascaded Filter (Red) and Individual Biquads (Green)');

grid on;

zoom on;

figure(2);

clf;

plot(f,h,'g--',f,h_cascaded,'r'); % plot frequency responses

axis([H_Min Nyquist_Freq V_Min V_Max]); % set the plotting range

xlabel('Frequency (Hz)');

ylabel('Magnitude (dB)');

title('Frequency Responses of Cascaded Filter (Red) and Individual Biquads (Green)');

grid on;

zoom on;



% Now write the fixed-point biquad filter coefficients to output *.h file

% Output format selection

hex_output = 1; % 0, 1

%eq_index = '_tx';

eq_index = '_rx';

if SF == 16000,

    eq_index = [eq_index,'_16kHz'];

elseif SF == 8000,

    eq_index = [eq_index,'_8kHz'];

end

filename = ['EQ_biquad_parameters',eq_index];

filename_n_sub = [filename,'.h'];

upper_filename = upper(filename);

file_header = ['__',upper_filename,'_H__'];

FID = fopen(filename_n_sub,'w');

fprintf(FID,'#ifndef %s\n',file_header);

fprintf(FID,'#define %s\n',file_header);

fprintf(FID,'/* This file is generated from EQ.m */\n\n');

fprintf(FID,'/* Sample Rate: %d */\n', SF);

fprintf(FID,'/* number of biquad sections: %d */\n\n', N_Stages);

fprintf(FID,'/* 32-bit EQ biquad filter coefficients in Q30 format */\n');

fprintf(FID,'/* numerator coefficients followed by 2nd and 3rd denominator coefficients */\n');

fprintf(FID,'/* one row per biquad section */\n');

%fprintf(FID,'INT32 coef_sos%s[%i] = {\n',eq_index, 5*N_Stages);

fprintf(FID,'#define EQ_coef_sos%s  \\\n',eq_index);

if hex_output == 0,

    for k=1:N_Stages-1,

        fprintf(FID,'    %11i, %11i, %11i, %11i, %11i,\\\n',coef_sos(k,:));

    end;

    fprintf(FID,'    %11i, %11i, %11i, %11i, %11i\n',coef_sos(N_Stages,:));

else

    for k=1:N_Stages,

        for j=1:5,

            if coef_sos(k,j) >= 0,

                fprintf(FID,'    0x%08x',coef_sos(k,j));

            else

                fprintf(FID,'    0x%08x',(double(coef_sos(k,j))+(2^(32))));

            end;

            if not((k==N_Stages) & (j==5)),

                % do not print , at the last line.

                fprintf(FID,',');

            end;

        end;

        if not((k==N_Stages) & (j==5)),

            % do not print \ at the last line.

            fprintf(FID,'\\\n');

        end;

    end;

end;



fprintf(FID,'\n\n');

fprintf(FID,'/* Per-biquad-section scaling gain coefficients in Q15 format */\n');

fprintf(FID,'/* one row per biquad section */\n');

fprintf(FID,'#define EQ_coef_g%s  \\\n',eq_index);

if hex_output == 0,

    for k=1:N_Stages-1,

        fprintf(FID,'    %6i, \\\n',coef_g(k));

    end;

        fprintf(FID,'    %6i\n',coef_g(N_Stages));

else

    for k=1:N_Stages,

        if coef_g(k) >= 0,

            fprintf(FID,'    0x%04x',coef_g(k));

        else

            fprintf(FID,'    0x%04x',(coef_g(k)+(2^16)));

        end;

        if not(k == N_Stages),

            % do not print \ at the last line.

            fprintf(FID,',\\\n');

        end;

    end;

end;

fprintf(FID,'\n\n');

fprintf(FID,'/* Gain shift amount applied after the last section. Range: -48..15 */\n');

if hex_output == 0,

    fprintf(FID,'#define EQ_gain%s  %i\n\n',eq_index, gain);

else

    if coef_g(k) >= 0,

        fprintf(FID,'#define EQ_gain%s    0x%04x\n\n',eq_index, gain);

    else

        fprintf(FID,'#define EQ_gain%s    0x%04x\n\n',eq_index, (gain+(2^16)));

    end;

end;



fprintf(FID,'#endif /* %s */\n',file_header);



fclose(FID);