|
|
function oF = design_biquad (iFs, iType, iFrequency, iBandGain, iThirdParam);
% This function designs the biquad filter, basing on filter type, sampling
% frequency, band gain, and Q-factor
frequency = iFrequency;
bandgain = iBandGain;
thirdparam = iThirdParam;
if frequency < 1.0,
frequency = 1.0;
end;
if frequency > 0.99*iFs/2.0
frequency = 0.99*iFs/2.0;
end;
omega = 2.0*pi*frequency/iFs;
cs = cos(omega);
sn = sin(omega);
A = 10.0^(bandgain/40.0);
alpha = sn/(2.0*thirdparam); % thirdparam = Q
switch upper(iType)
case 'PEAK' % peakingEQ
b0 = 1.0 + alpha*A;
b1 = -2.0*cs;
b2 = 1.0 - alpha*A;
a0 = 1.0 + alpha/A;
a1 = -2.0*cs;
a2 = 1.0 - alpha/A;
case 'LOW_SHELF' % lowShelf
b0 = A*( (A + 1.0) - (A - 1.0)*cs + 2.0*(A^0.5)*alpha );
b1 = 2.0*A*( (A - 1.0) - (A + 1.0)*cs );
b2 = A*( (A + 1.0) - (A - 1.0)*cs - 2.0*(A^0.5)*alpha );
a0 = (A + 1.0) + (A - 1.0)*cs + 2.0*(A^0.5)*alpha;
a1 = -2.0*( (A - 1.0) + (A + 1.0)*cs );
a2 = (A + 1.0) + (A - 1.0)*cs - 2.0*(A^0.5)*alpha;
case 'HIGH_SHELF' % highShelf
b0 = A*( (A + 1.0) + (A - 1.0)*cs + 2.0*(A^0.5)*alpha );
b1 = -2.0*A*( (A - 1.0) + (A + 1.0)*cs );
b2 = A*( (A + 1.0) + (A - 1.0)*cs - 2.0*(A^0.5)*alpha );
a0 = (A + 1.0) - (A - 1.0)*cs + 2.0*(A^0.5)*alpha;
a1 = 2*( (A - 1.0) - (A + 1.0)*cs );
a2 = (A + 1.0) - (A - 1.0)*cs - 2.0*(A^0.5)*alpha;
case 'LPF2' % second order LPF
b0 = (1.0 - cs)*0.5;
b1 = (1.0 - cs);
b2 = (1.0 - cs)*0.5;
a0 = 1.0 + alpha;
a1 = -2.0*cs;
a2 = 1.0 - alpha;
case 'HPF2' % second order HPF
b0 = (1.0 + cs)*0.5;
b1 =-(1.0 + cs);
b2 = (1.0 + cs)*0.5;
a0 = 1.0 + alpha;
a1 = -2.0*cs;
a2 = 1.0 - alpha;
case 'BPF' % BPF with constant 0 dB peak gain (for BPF with constant skirt gain and peak gain = Q: alpha->sn/2)
b0 = alpha;
b1 = 0.0;
b2 = -alpha;
a0 = 1.0 + alpha;
a1 = -2.0*cs;
a2 = 1.0 - alpha;
case 'NOTCH' % notch
b0 = 1.0;
b1 = -2.0*cs;
b2 = 1.0;
a0 = 1.0 + alpha;
a1 = -2.0*cs;
a2 = 1.0 - alpha;
case 'APF2' % APF second order
b0 = 1.0 - alpha;
b1 = -2.0*cs;
b2 = 1.0 + alpha;
a0 = 1.0 + alpha;
a1 = -2.0*cs;
a2 = 1.0 - alpha;
case 'APF1' % APF first order
% tbd. verify correct operation (check phase response), actually
% only listening test with sum of phased and unphased signal done,
% seems o.k.
omega= pi - omega; % tbd. verify this try&error hack
b0 = -1.0/( cos(omega) + sin(omega)/tan(pi/4.0 - omega/2.0) );
b1 = -1.0;
b2 = 0.0;
a0 = 1.0;
a1 = -b0;
a2 = 0.0;
case 'LPF1' % first order LP
b0 = 1.0/(1.0 + cos(pi*frequency/iFs)/sin(pi*frequency/iFs));
b1 = b0;
b2 = 0.0;
a0 = 1.0; % needed as later eliminated
a1 = 2.0*b0 - 1.0;
a2 = 0.0;
case 'HPF1' % first order HP
b0 = (cos(pi*frequency/iFs)/sin(pi*frequency/iFs))/(1.0 + cos(pi*frequency/iFs)/sin(pi*frequency/iFs));
b1 = -b0;
b2 = 0.0;
a0 = 1.0; % needed as later eliminated
a1 = 1.0 - 2.0*b0;
a2 = 0.0;
case 'PASS' % Pass-thru
b0 = 1.0;
b1 = 0.0;
b2 = 0.0;
a0 = 1.0; % needed as later eliminated
a1 = 0.0;
a2 = 0.0;
otherwise
error('VK: design_biquad internal error> indefinite type of filter');
end;
oF.A = [a0 a1 a2]/a0;
oF.Type = upper(iType);
oF.Fs = iFs;
oF.Frequency = iFrequency;
oF.ThirdParam = iThirdParam;
switch upper(iType)
case {'NOTCH', 'PEAK', 'HIGH_SHELF', 'LOW_SHELF'}
oF.B = [b0 b1 b2]/a0;
oF.BandGain_dB = 0;
oF.BandGainLin = 1;
otherwise
oF.B = A*[b0 b1 b2]/a0;
oF.BandGainLin = A;
oF.BandGain_dB = iBandGain;
end;
*******************************************************************
Пример использования
% This program demonstrates fading of one biquad in terms of frequency
% response change.
clear;
% c_Fs - sampling frequency, Hz
c_Fs = 44100;
% c_Num_Cascades - number of cascades of biquads
% c_Num_Cascades = 3;
% set the initial filter
% Fs, type, frequency, bandgain, thirdparam
% F0{1} = design_biquad(c_Fs, 'LPF2', 500, 0, 0.707);
% % set the target filter
% FN{1} = design_biquad(c_Fs, 'HPF2', 600, 0, 0.707);
%
% F0{2} = design_biquad(c_Fs, 'NOTCH', 200, 0, 6);
% FN{2} = design_biquad(c_Fs, 'PEAK', 1000, 6, 6);
% F0{1} = design_biquad(c_Fs, 'HPF2', 700, 0, 3);
% % set the target filter
% FN{1} = design_biquad(c_Fs, 'HPF1', 900, 0, 0.707);
% % FN{1} = design_biquad(c_Fs, 'LPF2', 7000, 0, 0.707);
%
% F0{2} = design_biquad(c_Fs, 'NOTCH', 2700, 0, 2);
% FN{2} = design_biquad(c_Fs, 'PEAK', 2760, 12, 2);
%
% F0{3} = design_biquad(c_Fs, 'PEAK', 6800, 6, 6);
% FN{3} = design_biquad(c_Fs, 'PASS', 4950, 0, 6);
%
% F0{4} = design_biquad(c_Fs, 'LPF1', 8000, 0, 0.707);
% FN{4} = design_biquad(c_Fs, 'APF2', 8000, 0, 0.707);
% For figure 3a-c
F0{1} = design_biquad(c_Fs, 'PEAK', 1000, 6, 3);
% set the target filter
FN{1} = design_biquad(c_Fs, 'LPF2', 7000, 0, 0.707);
% % For figure 2a-c
% F0{1} = design_biquad(c_Fs, 'HPF2', 100, 0, 0.707);
% % set the target filter
% FN{1} = design_biquad(c_Fs, 'LPF2', 1000, 0, 0.707);
% c_Num_Cascades - number of cascades of biquads
c_Num_Cascades = length(FN);
% c_mode_den - mode of denominator fading:
% =0 - linear fading;
% =1 - fading of poles
Params.c_mode_den = 1;
% c_mode_num - mode of numerator fading:
% =0 - fading via roots
% =1 - linear fading
% =2 - power fading
% =3 - parametric fading using normalization criteria ("sliding edges")
Params.c_mode_num = 3;
% c_eps - calculation accuracy
Params.c_eps = 0.00001;
% c_N - number of fading steps
Params.c_N = 10;
% c_FREQ_LEN - length of the frequency response to be plotted
c_FREQ_LEN = 4410;
% c_freq_units - type of the frequency units at the plot
c_freq_units = 'Hz';
% c_Dirac - Dirac's vector of length c_FREQ_LEN
c_Dirac = [1 zeros(1, c_FREQ_LEN-1)];
%col = cell(c_COLN, 1);
col_col = {'k'; 'b'; 'g'; 'c'; 'm'; 'r'; 'k:'; 'b:'; 'g:'; 'c:'; 'm:'; 'r:'};
col_bw = {'k', 'k:', 'k-.', 'k--'};
col = col_col;
% c_COLN - number of colors
c_COLN = length(col);
for I = 1: c_Num_Cascades
[Fi{I}.Bi, Fi{I}.Ai, Fi{I}.Gi] = fbqd_fading(F0{I}, FN{I}, Params);
end;
% Ai
% Bi
figure(2);
hold;
switch c_freq_units,
case 'Hz'
freq_vec = [0:floor(c_FREQ_LEN/2)-1]*c_Fs/c_FREQ_LEN;
case 'normalized_2pi'
freq_vec = [0:floor(c_FREQ_LEN/2)-1]*2*pi/c_FREQ_LEN;
case 'normalized'
freq_vec = [0:floor(c_FREQ_LEN/2)-1]/c_FREQ_LEN;
otherwise
error('Main> unsupported type of c_freq_units');
end;
for (I = 1:1:Params.c_N+1)
h = c_Dirac;
for (J = 1: c_Num_Cascades)
h = filter(Fi{J}.Bi(I, :), Fi{J}.Ai(I, :), h);
end;
H = abs(fft(h));
plot(freq_vec, H(1:floor(c_FREQ_LEN/2)), col{(mod(I-1, c_COLN) + 1)} );
end;
grid;
|
|
