Least P-norm optimal IIR Arbitrary magnitude filter design

I am designing multiple arbitrary magnitude response filter using iirlpnorm function. All the required magnitude response for these filters follow a pattern as they are all derived from same fitting function. But somehow the respone of the designed filters does not follow the original trend of the required response. I am also liitle confused with how the weights for each frequency points effect the filter response as slight variation in weights is seen to severely distort the response of the filter.Attached is the response required.All the response are obtained from a power fit. Any help or clarification regarding iirlpnorm function will would be highly appreciated.

Thanks

[Mathieu NOE]

do you have a working code to share which shows what you want to achieve ?

[rizal_m]

Here is the code. The data pertains to the green line (one location) in the graph attached. I can design the filter for each individual plots in the graph (represents filter property at different location) but they don't follow the trend as seen in the graph of required response.

Fs = 20000000;  % Sampling Frequency

Nb   = 1;                                        % Numerator Order
Na   = 6;                                        % Denominator Order
E    = [0 475000];                               % Frequency Edges
W    = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];  % Weight Vector
P    = [2 128];                                  % P'th norm
dens = 20;                                       % Density Factor

% Frequency Vector
F    = [0 50000 75000 100000 125000 150000 175000 200000 225000 250000 275000 300000 325000 350000 375000 400000 425000 450000 475000];

% Amplitude Vector
A    = [1  0.418142 0.2421 0.134898 0.072923 0.038437 0.01982 0.010028 0.004986 0.002440 0.001177 0.000560 0.000263 0.000122 5.625e-05 2.55e-05 1.151e-05 5.13e-06 2.27e-06];

[b,a,err,sos_var,g] = iirlpnorm(Nb, Na, F/(Fs/2), E/(Fs/2), A, W, P, ...
                                {dens});
Hd                  = dfilt.df2sos(sos_var, g);

[Mathieu NOE]

hello again

you mention a magnitude vs frequency , but don't you have also a phase associated with it ?

[rizal_m]

Hello,

I dont have a phase response associated with the data I am working with. Only the amplitude attenuation and the response patterns as seen in the graph are of interest. I hope this helps.

Thanks

[Mathieu NOE]

may I ask you where does this magnitude data come from ? how was it generated ?

[rizal_m]

it was obatined from experiments. The magnitude data represents attenuation due to damping.

[Mathieu NOE]

so you measured some kind of transfer function ? then you should also have the phase information somehow

[rizal_m]

hello,

like I mentioned above the we are not interested in the phase response. We are just measuring amplitude drop of different central frequencies irrespective of the phase.I just want the filter amplitude response to match the experimental data.

[Mathieu NOE]

hello

i am only half way ... without the help of the DSP Tbx taht I don't have

so I am working with other means

I coud get a FIR model working but I am not yet ready to show a IIR solution...

Fs = 20000000;  % Sampling Frequency
Nb   = 1;                                        % Numerator Order
Na   = 6;                                        % Denominator Order
E    = [0 475000];                               % Frequency Edges
W    = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];  % Weight Vector
P    = [2 128];                                  % P'th norm
dens = 20;                                       % Density Factor
% Frequency Vector
F    = [0 50000 75000 100000 125000 150000 175000 200000 225000 250000 275000 300000 325000 350000 375000 400000 425000 450000 475000];
% Amplitude Vector
A    = [1  0.418142 0.2421 0.134898 0.072923 0.038437 0.01982 0.010028 0.004986 0.002440 0.001177 0.000560 0.000263 0.000122 5.625e-05 2.55e-05 1.151e-05 5.13e-06 2.27e-06];

% interpolate F,A data 
Ni = 300;
Fi = linspace(min(F),max(F),Ni);
Ai  = interp1(F,A,Fi);

% FIR approximation (magnitude match only)
N = 301;
FIR = fsamp( Ai, N ); % NB : filter length N must not be greater than 2*(length(mag)-1) = 2*Ni-1
H = freqz(FIR,1,N);
freq2 = linspace(min(F),max(F),N);

figure(1)
% subplot(211)
plot(F,20*log10(A),freq2,20*log10(abs(H))); % y scale in dB (20log10())
% subplot(212)
% plot(freq1,zeros(1,Ni),freq2,180/pi*(angle(H))); % y scale in degrees

% % IIR filter design
% NB = 1; % numerator order
% NA = 6; % denominator order
% % W = 2*pi*freqs/fs;
% W = pi*freq_normalized2;
% Wt = ones(size(W));
% % Wt = linspace(1,0.1,length(W));
% TOL = 1e-2;
% ITER = 1e4;
% [B,A] = invfreqz(H,W,NB,NA,Wt,ITER,TOL);
% hverif = freqz(B,A,N);
% 
% figure(2)
% subplot(211),plot(freq_normalized1,20*log10(Ai),freq_normalized2,20*log10(abs(H)),freq_normalized2,20*log(abs(hverif)));% y scale in dB (20log10())
% subplot(212),plot(freq_normalized1,zeros(1,Ni),freq_normalized2,180/pi*(angle(H)),freq_normalized2,180/pi*(angle(hverif))); % y scale in degrees
% legend('spec','FIR realisation','IIR realisation'); 
% xlabel('Frequency (Hz)');
% ylabel('modulus (dB)');
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% help found here : 
% https://dsp.stackexchange.com/questions/87007/how-to-create-a-filter-that-matches-a-desired-magnitude-response
function h = fsamp(mag,N)
% function h = fsamp(mag,N)
% frequency sampling design of linear phase FIR filter
%
% mag ... desired magnitude response on an equidistant grid
%         between 0 and Nyquist
% N   ... filter length (N <= 2*length(mag) - 1 )
% 
% for even N the frequency grid is defined by (M = length(mag)):
%
% f_k = k/(M-1), k=0,1,...,M-1 (includes Nyquist = 1)
%
% for odd N the grid is f_k = 2/(2*M-1), k=0,1,...,M-1 (Nyquist is not included!)
%
% M. Lang, [e-mail]

M = length(mag);
Neven = (rem(N,2) == 0);
if Neven,
    n = 2*(M-1);                    % FFT length
    if ( N > n ),
        error('fsamp: filter length N must not be greater than 2*(length(mag)-1).');
    end
    w = (0:M-1)'/(M-1)*pi;          % frequency grid
    D = mag(:).*exp(-1i*w*(n-1)/2); % desired complex frequency response
    D = [D;conj(D(M-1:-1:2))];      % extend to 2*Nyquist
else
    n = 2*M-1;
    if ( N > n ),
        error('fsamp: filter length N must not be greater than 2*(length(mag)-1).');
    end
    w = (0:M-1)'*2*pi/(2*M-1);
    D = mag(:).*exp(-1i*w*(n-1)/2);
    D = [D;conj(D(M:-1:2))];
end

% impulse response via IFFT
h = real(ifft(D));

% windowing (Hamming)
I = (n-N)/2+1;
ww = hamming(N);
h = h(I:I+N-1).*ww(:);

end

[Mathieu NOE]

this is an attempt to optimize a IIR filter to your data

with a order = 6 it's not convincing, I can get fairly good results above order 16

now I cannot compare with yur results as I don't have the dsp tbx

result

code

Fs = 20000000;  % Sampling Frequency
Nb   = 1;                                        % Numerator Order
Na   = 6;                                        % Denominator Order
E    = [0 475000];                               % Frequency Edges
W    = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];  % Weight Vector
P    = [2 128];                                  % P'th norm
dens = 20;                                       % Density Factor
% Frequency Vector
F    = [0 50000 75000 100000 125000 150000 175000 200000 225000 250000 275000 300000 325000 350000 375000 400000 425000 450000 475000];
% Amplitude Vector
A    = [1  0.418142 0.2421 0.134898 0.072923 0.038437 0.01982 0.010028 0.004986 0.002440 0.001177 0.000560 0.000263 0.000122 5.625e-05 2.55e-05 1.151e-05 5.13e-06 2.27e-06];

% interpolate F,A data 
Ni = 100;
Fi = linspace(min(F),max(F),Ni);
AdB  = interp1(F,20*log10(A),Fi);

max_order = 16;

IC = linspace(max(Fi)/max_order,max(Fi),max_order);
po = -2*pi*IC;
ze = zeros(numel(po),1);

[~,a] = zp2tf(ze,po,1);
b= zeros(size(a));
be = prod(po);
b(end) = be;

H = freqs(b,a,2*pi*Fi);
H_dB_initial = 20*log10(abs(H));

%% optimisation fminsearch
global AdB Fi
x = fminsearch(@objectivefcn1,IC);
x = abs(x); % force x to be positive
po = -2*pi*x; % poles (must be negative)
ze = zeros(numel(po),1); %zeroes init

[~,a] = zp2tf(ze,po,1);
b= zeros(size(a));
b(end) = prod(po); % so TF static gain = 1

H = freqs(b,a,2*pi*Fi);
H_dB_final = 20*log10(abs(H));

figure(1)
plot(Fi,AdB,Fi,H_dB_final); % y scale in dB (20log10())
xlabel('Freq (Hz)');
ylabel('Gain (dB)');

legend('specification','realisation');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = objectivefcn1(x,AdB,Fi)
global AdB Fi
    x = abs(x); % force x to be positive
    po = -2*pi*x;
    ze = zeros(numel(po),1);

    [~,a] = zp2tf(ze,po,1);
    b= zeros(size(a));
    be = prod(po);
    b(end) = be;

    H = freqs(b,a,2*pi*Fi);
    H_dB = 20*log10(abs(H));

    f = norm(H_dB - AdB);

end

[rizal_m]

Hello,

Thank you for the code. It works perfectly for higher frequency components and the response also matches the required response. But it does not behave well in lower frequencies. Like I stated above, each set of data represents drop in amplitude at different distance. So we would expect the amplitude drop of shorter distance to be less than longer one. In other words, the response of longer distnace when plotted together with shorter distance should be stacked just below shorter distance. But they seem to cross each other at lower frequencies as seen in Figure below. Is there any way where we can reduce the error or optimize the filter even more at lower frequencies?

Thanks

[Mathieu NOE]

hello again

I made a few attempts to get better results at lower frequencies

In order to have a different view on the results, I made a second plot with x in log scale (as usually we do Bode plots with log axis both in x and y)

also you noticed that the optimization is done on interpolated data, not the raw data (to improve results) . Notice also that I did a linear interpolation on dB values (and not vive versa, which would have been taking the 20*log10 of A interpolated values (linear) , because this option would lead to straigth segments in the Bode plot and that is not trully representative of the behaviour of a filter).

now , to have more points in the lower frequency range, I switched from linearly spaced x values to log spaced values ; this can slightly improve things

also I opted for order = 21 which allows a better low frequency match

as you will notice on the second plot (with log spaced x values) , we lack some input data as the only first frequency is already at a fairly high value (50000 Hz) , so can we have more points between 0 and 50000 Hz ?

last point : to avoid log of 0 , I changed the first frequency value to a non zero value (here 1, could be also 10 or 100 as this remains quite low (quasi static) if we compare to the next value = 50000 Hz

So , this is what we get with order = 21 , second plot (log scale on x axis) , and data are linearly interpolated (vs freq values F) : Fi = linspace(min(F),max(F),Ni);

And this is what we get with order = 21 , second plot (log scale on x axis) , and data are log interpolated (vs freq values F) : Fi = logspace(log10(min(F)),log10(max(F)),Ni);

it's better !!

again , do you have intermediate values between F = 0 and 50000 Hz , so we can better cover this range ?

[rizal_m]

hello,

We don't have any data on that range and in fact most of the data was extrapolated from the data we had from the experiments.

Thanks

[Mathieu NOE]

ok

does at least the proposed code make a better job in terms of avoiding curves to cross each others ?

[rizal_m]

yeah, I think it might work. I might have to tweak the interpolated frequency individually for each distance. I am working on it. Thank you for the help.

[Mathieu NOE]

ok let me know how things turn

[Mathieu NOE]

so far so good?

[rizal_m]

Hi,

Thanks for the follow-up. I think it should be good. I tried both the FIR and irr p-norm. The response is satisfactory in the freuqnecy band we are interested in.

Thanks for the help.