How to get magnitude values for a Parametric EQ?

I was doing pretty well with analog filters, but I can't seem to get similar results with digital filters.
Here's what I tried first. Link to my xlsx data.
%% load stuff
M = readtable('dB Technologies T4 PS1.xlsx','Sheet','TF');
%% create filter
Fs = 96e3;
N1 = 2;
G = -5; % 5 dB
Wo = 100/(Fs/2);
BW = 4000/(Fs/2); % Need to figure out how to insert BW in octaves.
[b,a] = designParamEQ(N1,G,Wo,BW);
%% calculate magnitude response
w = M.Frequency_Hz;
h = freqz(b,a,w);
magdB = mag2db(abs(h));
That didn't work, so I tried create a biquad filter from the parametric filter.
BQ = dsp.BiquadFilter('SOSMatrix',[b.',[1,a.']]);
h = freqz(BQ);
That didn't work either. I thought maybe it had something to do with my frequency needing to be converted to radians per sample...
wRad = w*((2*pi)/(Fs/2));
But that didn't work either. What am I missing here?
By the way, I was able to visualize the BQ filter with fvtool, so I could see it was working. I just have been able to get a magnitude response form it.

 Risposta accettata

Brian Hemmat
Brian Hemmat il 22 Ago 2020
Modificato: Brian Hemmat il 27 Ago 2020
Hi Nathan,
The syntax of freqz that you were using in your code above expects w to be normalized fequency. I ran your code and its not normalized. You can specify frequency in Hz by also specifying the sample rate.
Another issue (updated answer): is that the default orientation of the coefficients output from designParamEQ was not the orientation expected by fvtool and freqz.
M = readtable('spec.xlsx','Sheet','TF');
Fs = 96e3;
N = 2;
G = -5; % 5 dB
Wo = 1000/(Fs/2); % center frequency
BW = 1; %bandwidth in octaves
Q = sqrt(2^BW)/((2^BW)-1); % convert bandwidth to Q
BWnorm = Wo/Q; % then to normalized bandwidth
[b,a] = designParamEQ(N1,G,Wo,BWnorm,'Orientation','row');
fvtool([b,a])
%%
% Magnitude and phase of parametric EQ
f = M.Frequency_Hz;
h = freqz([b,a],f,Fs);
magdB = mag2db(abs(h));
% plot
semilogx(f,magdB) % the filter
grid on
% ylim([-24 24]);
% xlim([20 20000]);
xlabel('Frequency')
ylabel('Magnitude')
legend('Mag Parametric EQ','Location','best');

4 Commenti

Brian, thank you so much for this. I marked the question as answer because I think you did, at 10kHz, but if you lower the center frequency to 1kHz, it falls apart. It seems like maybe bandwidth does not work like I expect it to work.
This code, for example.
%% create Parametric filter
Fs = 96e3;
N1 = 2;
G = -5; % 5 dB
Wo = 1000/(Fs/2); % center frequency
BW = 1; %bandwidth in octaves
Q = sqrt(2^BW)/((2^BW)-1); % convert bandwidth to Q
BWnorm = Wo/Q; % then to normalized bandwidth
[b,a] = designParamEQ(N1,G,Wo,BWnorm);
% Magnitude and phase of parametric EQ
f = M.Frequency_Hz;
h = freqz(b,a,f,Fs);
magdB = mag2db(abs(h));
% plot
subplot(2,1,1)
semilogx(f,magdB) % the filter
grid on
ylim([-24 24]);
xlim([20 20000]);
xlabel('Frequency')
ylabel('Magnitude')
legend('Mag Parametric EQ','Location','best');
Creates this plot. I'm expecting to see a nice 5dB dip at 1kHz.
I also noticed in your plots that the filter is at -25dB, not -5dB.
Thanks for any more guidance you might have.
Hey Brian, it would be great to get more feedback on this, or maybe a resource where I could learn more about what's happening.
Hi Nathan,
Sorry, I led you astray on my first answer. There was a second issue that the orientation of the coefficients output from the design function was not compatible with fvtool and freqz. Consult the documentation for more details:
We did it! Well, you did it. :)
Brian, thank you so much. It's been a struggle for me to get started because I know a little about being a sound engineer, but almost nothing about math or programming.
Here's the result.
Would you be willing to take a stab at helping figure out the same problem, but with an allpass filter? https://www.mathworks.com/matlabcentral/answers/581223-why-is-my-1khz-allpass-filter-returning-1-125khz

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Measurements and Spatial Audio in Centro assistenza e File Exchange

Prodotti

Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by