Proper way to call designfilt for filtfilt

The command designfit is used to design a filter with user high level specification of the response funtion. The result of designfilter can be applied on a signal by using filter function.

In my caseI want to apply with zero-phase filter filtfilt instead of filter. Howeever this one squares the amplitude response function the the high level specification is not right.

Is there a proper way of compensate for this squaring in a high level? or any existing filter design function suitable for filtfilt?

Right now I use freqz function and adjust manually the parameter until te response curve looks empirically "right". I could do some more precise numerical optimisation to.

Do you kwon a better way?

[Paul]

Hi Bruno,

Can you provide an example call to designfilt for a particular case of interest?

There may be different options to explore depending on if the filter is IIR or FIR. In particular, if FIR then lowpass or highpass, etc. maybe of interest. IIRC, those functions "compensate for the delay" without doing the forward/backward processing that changes the amplitude response.

[Bruno Luong]

Here is a simple example

CutoffFrequency = 1500;
facq = 2e5;

Vfilter = designfilt('lowpassfir', 'FilterOrder', ceil(facq/CutoffFrequency), ...
    'CutoffFrequency', CutoffFrequency, 'SampleRate', facq);

% This responnse corresponds to filter not filtfilt
% I want to "match" filtfilt response to this (especially in the zone 0-4kHz, or magnitude > -50 dB)
fig = figure('Name', 'Filter caracteristic');
ax3 = axes('Parent', fig);
[h, f] = freqz(Vfilter, 1024, facq);
semilogx(ax3, 1e-3*f, 20*log10(abs(h)/max(h)))
xlabel(ax3, 'f [kHz]');
ylabel(ax3, 'Magnitude [dB]');
grid(ax3, 'on');
title(ax3, 'Filter magnitude response');

[Star Strider]

If you already know the exact transfer function you want, experiment with firls. (In designfilt, use 'DesignMethod','ls' to do that, along with 'Frequencies' and 'Amplitudes' vectors, and a 'Weights' vector if you want to assign them.) It can very closely approximate most desired transfer functions, and it allows a weighting vector to give more weight to specific frequencies. If you have a long signal, you can use a higher-order filter (increased precision, decreased computational efficiency).

I rarely use designfilt, so I have llittle experience with it, preferring command-line filter design.

[Mathieu NOE]

hello

why not simply reduce the filter order by a factor 2 when you want to use filtfilt instead of filter ? of course that makes sense only if the filter command is used for even order filters

[Bruno Luong]

@Star Strider I like the idea, it seems a priori more generic than Paul method, but actually it also applied for FIR and more complicated.

Answer 1

[Paul]

The example above is a low pass, linear phase FIR filter. Maybe a simple delay as is implemented in lowpass is acceptable?

CutoffFrequency = 1500;
facq = 2e5;
Vfilter = designfilt('lowpassfir', 'FilterOrder', ceil(facq/CutoffFrequency), ...
    'CutoffFrequency', CutoffFrequency, 'SampleRate', facq);

filter a sine wave

f0 = 2e3;
n = 0:1000;
t = n/facq;
x = sin(2*pi*f0*t);
y1 = filter(Vfilter.Coefficients,1,x);
figure
plot(t,x,t,y1),grid

The amplitude and phase shift are apparent consistent with the frequency response of Vfilter.

Because the linear phase FIR filter has a constant delay, we can zero pad the end of the input by the number of delayed samples, then filter, then delete number of delayed samples from the front of the output.

ndelay = Vfilter.FilterOrder/2;
y2 = filter(Vfilter.Coefficients,1,[x zeros(1,ndelay)]);
y2(1:ndelay) = [];
figure
plot(t,x,t,y2),grid

Same steady "state amplitude" of the output, but no more steady state phase delay. Here "steady state" is referring to the "middle of the response" as the effect of the impulse response of the filter is evident on the first and last cycles of the output.

Use a frequency sweep to show the zero-phase over all frequencies.

x = sweeptone(6,2,facq,SweepFrequencyRange=[100 facq/2]).';
y3 = filter(Vfilter.Coefficients,1,[x zeros(1,ndelay)]);
y3(1:ndelay) = [];

[H,fvec] = tfestimate(x,y3,[],[],[],facq);
figure
subplot(211)
semilogx(1e-3*fvec,db(abs(H)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
xlim([0.1 100])
subplot(212);
semilogx(1e-3*fvec,180/pi*angle(H),1e-3*fvec,180/pi*angle(freqz(Vfilter,fvec,facq)));grid
xlim([0.1 100])

I'm pretty sure that the blue phase starts getting hairy because of rounding errors accumulating with 134-tap filter and then trying to take the angle of very small complex numbers.

[Bruno Luong]

Thanks. The technique seems similar to this https://www.mathworks.com/help/signal/ug/compensate-for-the-delay-introduced-by-an-fir-filter.html

Also seems joints @Mathieu NOE idea of reduce the filter-order by half (I'm not sure yet however the correct way to implement it).

[Bruno Luong]

So it seems beyond FIR there is no straighforward method?

[Paul]

Yes, essentially the same as the linked doc page on compensating for the delay introduced by the FIR filter, at least for a filter of even order as in this case.

I'll give some more thought to the IIR filter case. Can you give an example of an IIR filter to work with, preferably with a call to designfilt?

Code below is essentially the same as above, but computes the relevant frequency responses without discarding the elements in the output that are shifted into negative time (and I used a simple FFT ratio to compute the transfer function frequency response instead of tfestimate). Keeping those points makes the frequency response cleaner (even w/o the the averaging that tfestimate does).

CutoffFrequency = 1500;
facq = 2e5;
Vfilter = designfilt('lowpassfir', 'FilterOrder', ceil(facq/CutoffFrequency), ...
    'CutoffFrequency', CutoffFrequency, 'SampleRate', facq);

Nh = Vfilter.FilterOrder;
Ndelay = Nh/2;

x = sweeptone(6,2,facq,SweepFrequencyRange=[100 facq/2]).';
Nx = numel(x); % length of input

y3 = filter(Vfilter.Coefficients,1,[x zeros(1,2*Nh)]);
ny3 = (0:numel(y3)) - Ndelay;  % time of output

% DTFT of the input
[X,w]  = freqz(x,1,4096);
% DTFT of the output
Y3     = freqz(y3,1,w).*exp(-1j*ny3(1)*w);

H = Y3./X; 
fvec = w*facq/2/pi;

figure
subplot(211)
semilogx(1e-3*fvec,db(abs(H)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
xlim([0.1 100])
subplot(212);
semilogx(1e-3*fvec,180/pi*angle(H),1e-3*fvec,180/pi*angle(freqz(Vfilter,fvec,facq)));grid
xlim([0.1 100])

[Bruno Luong]

Here is an example of IIR. I use typically either this of the earlier FIR in my application:

PassbandFrequency = 1000;
FilterOrder = 5;
StopbandAttenuation = 60;
facq = 2e5;
Vfilter = designfilt('lowpassiir', 'FilterOrder', FilterOrder, ...
    'PassbandFrequency', PassbandFrequency, ...
    'StopbandAttenuation', StopbandAttenuation, ...
    'SampleRate', facq);

% This responnse corresponds to filter not filtfilt
% I want to "match" filtfilt response to this (especially in the zone 0-4kHz, or magnitude > -50 dB)
fig = figure('Name', 'Filter caracteristic');
ax3 = axes('Parent', fig);
[h, f] = freqz(Vfilter, 1024, facq);
semilogx(ax3, 1e-3*f, 20*log10(abs(h)/max(h)))
xlabel(ax3, 'f [kHz]');
ylabel(ax3, 'Magnitude [dB]');
grid(ax3, 'on');
title(ax3, 'IIR Filter magnitude response');

[Paul]

I tried using designfilt with the arbmagfir response to attempt to match the sqrt of the amplitude response of VFilter, but didn't get very far (and, frankly, I didn't try very hard).

I then switched to a method base on the forward-backward conv shown in this thread and got pretty close, but not close enough to leave here. I'll come back to this when I have a chance, but it might be a while.

Did you try asking over on dsp.stackexchange.com?

[Bruno Luong]

@Paul, Thanks. Don't worry I do manually and it's OK for now. I though the topic is intesting and challenging. So come back whenever you have something new to share.

I'm not subscibed to stackexchange.

[Paul]

Hi Bruno,

I'll leave this here for now.

The basic idea is to estimate the impulse response of the IIR filter that has magnitude equal to the sqrt of the design filter. Then, with that impulse response do forward-backward filtering in the time domain on the input signal. Frankly, I was disappointed the method doesn't work better for this example, but I didn't have time to explore any further. Here we go ....

Define the design filter. The goal is operate on an input signal such that transfer function to the output has the same magnitude as the design filter, but has zero phase.

PassbandFrequency   = 1000;
FilterOrder         = 5;
StopbandAttenuation = 60;
facq                = 2e5;
Vfilter = designfilt('lowpassiir', 'FilterOrder', FilterOrder, ...
    'PassbandFrequency', PassbandFrequency, ...
    'StopbandAttenuation', StopbandAttenuation, ...
    'SampleRate', facq);

Compute and plot the impulse response of the deisgn filter

Vimp = impz(Vfilter);
NVimp = numel(Vimp);
figure
plot(Vimp);

For all intents and purposes, we'll assume that Vimp[n] = 0 for n >= NVimp.

Compute one period of the DTFT of Vfilter. I was playing around with the number of frequency samples.

N = 4096;
[Vfreq,wfreq] = freqz(Vfilter,N*4,'whole');
Vfreq(end+1) = Vfreq(1); wfreq(end+1) = 2*pi;

The frequency response with same phase as and sqrt of the magnitude of Vfilter is

sqrtVfreq = sqrt(abs(Vfreq)).*exp(1j*angle(Vfreq));

Define an anonymous function to evaluate the frequency response at arbitrary frequencies. I'm using linear interpolation here, which might not be good enough. I also wrote a short function at the bottom of the post to compute sqrtVfreq based on a call to freqz(Vfilter, ....), but it didn't seem to help much. So I'll go with this for now.

sqrtVfreqz = @(w) interp1(wfreq,sqrtVfreq,w);

Compute the impulse response from sqrtVfreqz using the IDTFT integral. Should use integral, but that was taking too long to run. So I'll just use trapz. The assumptions are that the sqrt filter is causal, and its impulse response is real, and it can be adequately represented by the same number of elements as the effective length of the impulse response of the design filter.

for ii = 1:NVimp
    
    n = ii-1;
    sqrtVimp(ii) = trapz(wfreq,sqrtVfreqz(wfreq).*exp(1j*n*wfreq))/2/pi;
    %sqrtVimp(ii) = trapz(wfreq,func(Vfilter,wfreq).*exp(1j*n*wfreq))/2/pi;
    
end

Add sqrtVimp to the plot to compare to Vimp. Its not that different (at this scale) and the assumption of when it's effectively zero seems reasonable.

hold on;plot(real(sqrtVimp))

Verify that the imaginary part is numerical trash.

figure
plot(imag(sqrtVimp))

Retain only the real part.

sqrtVimp = real(sqrtVimp);

Could also have checked for causality using the IDFT integral, but I haven't done that. Actually, causality of sqrtVimp might not be important (in the sense that we can account for it being noncausal), as long as it decays to close to zero in finite (negative) time. Hmm, need to think about this some more.

Define the input test signal

x = sweeptone(6,2,facq,SweepFrequencyRange=[100 facq/2]).';
x = x(:).';  % ensure row vector

The rest of this code was previously used in filtfilt Does Not Perform Zero-Phase Filtering. Does It Matter?

Rename the impulse response of the sqrt filter for convenience and ensure it's a row vector.

h = sqrtVimp(:).';              % sqrt filter impulse resonse

Use time domain convolution to implement the zero phase filter.

y2 = conv(fliplr(h),conv(h,x)); % zero-phase, filtered output

Time domain indices of y2. Note that y2 extends to negative time, as must (I think) be the case for zero-phase filtering a signal that starts at n = 0.

Nh = numel(h); % length of impulse response
Nx = numel(x); % length of input
ny2 = -(Nh-1) : Nh + Nx - 2;  % time of output

DTFT of the input signal (half-period)

[X,w]  = freqz(x,1,N);

DTFT of the output signal (half-period)

Y2     = freqz(y2,1,w).*exp(-1j*ny2(1)*w);

Corresponding frequency vector

fvec = w/2/pi*facq;

Plots for comparison of the FFT ratio of Y2/X and the original design filter.

figure
subplot(211)
semilogx(1e-3*fvec,db(abs(Y2./X)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid

subplot(212);
semilogx(1e-3*fvec,180/pi*angle(Y2./X),1e-3*fvec,180/pi*angle(freqz(Vfilter,fvec,facq)));grid

The zero phase of Y2/X looks good. The gain of Y2/X looks kind of similar to the magnitude of the design filter, but I really thought (hoped?) it would be better.

Zoom in on the gain at low frequecy.

figure
semilogx(1e-3*fvec,db(abs(Y2./X)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
xlim([0.1 1])

It kind of has a similar shape, but might not match well enough.

I might come back to this in the future. I think the initial areas for exploration would be to use the function below to evaluate sqrtVfreq, use something better than trapz for the IDTFT, and check that causality assumption on sqrtVimp.

Or maybe there is one or more fundamental errors in the math or the code or both.

Don't know if this method can be made to work better, but it was interesting to think about.

Disclaimer: I only tested this method on this Vfilter and with only this test input.

function sqrtVfreq = func(Vfilter,w)

sqrtVfreq = freqz(Vfilter,[-1;w(:)]); % ensure the w input has at least two elements
sqrtVfreq = sqrt(abs(sqrtVfreq)).*exp(1j*angle(sqrtVfreq));
sqrtVfreq(1) = [];
sqrtVfreq = reshape(sqrtVfreq,size(w));

end

[Bruno Luong]

Hi Paul. Thanks.

It seems to me that we want to approximate the square-root of the IIR norm of fractional transfert function by another norm of the ractional transfert function. Unless the original zeros and poles of the IIR are even order, exact solution won't exist, so only approxilation is possible.

I'm have impression that optimization will perform better than what you get above. Not sure why such tool is not implemented for IIR, but only for FIR.

[Paul]

The assumption that sqrtVimp is causal is incorrect and makes a big difference. Code that follows is almost the same as above, but computes sqrtVImp for negative time.

PassbandFrequency   = 1000;
FilterOrder         = 5;
StopbandAttenuation = 60;
facq                = 2e5;
Vfilter = designfilt('lowpassiir', 'FilterOrder', FilterOrder, ...
    'PassbandFrequency', PassbandFrequency, ...
    'StopbandAttenuation', StopbandAttenuation, ...
    'SampleRate', facq);

Vimp = impz(Vfilter);
NVimp = numel(Vimp);
N = 4096;
[Vfreq,wfreq] = freqz(Vfilter,N*4,'whole');
Vfreq(end+1) = Vfreq(1); wfreq(end+1) = 2*pi;

sqrtVfreq = sqrt(abs(Vfreq)).*exp(1j*angle(Vfreq));
sqrtVfreqz = @(w) interp1(wfreq,sqrtVfreq,w);

I changed this loop to compute sqrtVimp for negative time as well. The outer loop is now the discrete-time variable and the Matlab array index is computed in the loop.

for n = -NVimp:NVimp
    
    ii = n + NVimp + 1;
    sqrtVimp(ii) = trapz(wfreq,sqrtVfreqz(wfreq).*exp(1j*n*wfreq))/2/pi;
    %sqrtVimp(ii) = trapz(wfreq,func(Vfilter,wfreq).*exp(1j*n*wfreq))/2/pi;
    
end

figure
plot(-NVimp:NVimp,real(sqrtVimp)),grid

We see that sqrtVimp extends into negative time, so noncausal.

figure
plot(-NVimp:NVimp,imag(sqrtVimp))

sqrtVimp = real(sqrtVimp);

x = sweeptone(6,2,facq,SweepFrequencyRange=[100 facq/2]).';
x = x(:).';  % ensure row vector

h = sqrtVimp(:).';              % sqrt filter impulse resonse

y2 = conv(fliplr(h),conv(h,x)); % zero-phase, filtered output
Nh = numel(h); % length of impulse response
Nx = numel(x); % length of input

edit: 8 Nov 2023, corrected expression ny2, but it so happened that ny2(1) was correct in the original expression, and ny2(1) is the only term used.

Here, I had to adjust the discrete-time indices of y2 to account for the fact that sqrtVimp(1) corresponds to n = -2*NVimp in the time domain.

%ny2 = -(Nh-1) : Nh + Nx - 2 - NVimp;  % time of output
ny2 = (0:numel(y2)-1) - 2*NVimp;  % time of output

[X,w]  = freqz(x,1,N);
Y2     = freqz(y2,1,w).*exp(-1j*ny2(1)*w);
fvec = w/2/pi*facq;

figure
subplot(211)
semilogx(1e-3*fvec,db(abs(Y2./X)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
subplot(212);
semilogx(1e-3*fvec,180/pi*angle(Y2./X),1e-3*fvec,180/pi*angle(freqz(Vfilter,fvec,facq)));grid

figure
semilogx(1e-3*fvec,db(abs(Y2./X)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
xlim([0.1 1])

Those plots look very good.

Same disclaimers as above.

Comment: This approach is not specific to IIR, so should work for FIR as well (I didn't try it), though the simple shift method for FIR discussed previously is easier.

function sqrtVfreq = func(Vfilter,w)
sqrtVfreq = freqz(Vfilter,[-1;w(:)]); % ensure the w input has at least two elements
sqrtVfreq = sqrt(abs(sqrtVfreq)).*exp(1j*angle(sqrtVfreq));
sqrtVfreq(1) = [];
sqrtVfreq = reshape(sqrtVfreq,size(w));
end

[Bruno Luong]

@Paul magnificient ! You got it, hat off.

[Bruno Luong]

@Paul I have a question, you select the original phase for sqrtVfreq

sqrtVfreq = sqrt(abs(Vfreq)).*exp(1j*angle(Vfreq));

Is it a requirement or can you for example do time contraction

sqrtVfreq = sqrt(abs(Vfreq)).*exp(1j*angle(Vfreq) / 2); % divided the phase by half, or one of the sqrt(Vfreq)

?

[Paul]

I was wondering the same thing too, but didn’t pursue it. Only on mobile for a while, so can’t really do too much. IIRC, there is a requirement that sqrtVimp is real, and I thought that was likely to happen if I didn’t mess with the phase, but there may well be other options.

[Paul]

Hi Bruno,

I will hypothesize that one (the only?) requirement is that the sqrtVfreq has to be conjugate symmetric around omega = pi (rad/sample). In the code below, I'm comparing multiplicative factors on the phase of 0.5 and 0. It would be interesting to modify the phase by some other function of omega while maintaining conjugate symmetry (I haven't tried it).

The code also uses a different input, x, to allow it to run a bit faster.

Speaking of which, I suspect the vast majority of the exectution time is spent in the loop that approximates the IDFT. Maybe vectorizing that loop would be better, and, if we're willing to assert that the IDFT is real (it should be) then it would probably be faster to compute only the real part of the IDFT from the outset.

Also, I had an error in the computation of the ny2 vector, though ny2(1), which is the value used in the code happened to be correct (also corrected in the comment above).

PassbandFrequency   = 1000;
FilterOrder         = 5;
StopbandAttenuation = 60;
facq                = 2e5;
Vfilter = designfilt('lowpassiir', 'FilterOrder', FilterOrder, ...
    'PassbandFrequency', PassbandFrequency, ...
    'StopbandAttenuation', StopbandAttenuation, ...
    'SampleRate', facq);

N = 4096;

%x = sweeptone(6,2,facq,SweepFrequencyRange=[100 facq/2]).';
rng(100);
x = rand(1000,1);

phase_factor = [1 , 0.5 , 0];


for ii = 1:numel(phase_factor)

    disp(['phase_factor = ',num2str(phase_factor(ii))])
    [y2{ii},ny2] = analyze_this(Vfilter,N,x,phase_factor(ii));

end

phase_factor = 1

phase_factor = 0.5

phase_factor = 0

Looking at the results above, it appears that phase_factor = 0.5 did not work quite as well as with 0 or 1, probably because the range -NVimp:NVimp isn't long enough for the IIR response of sqrtVImp to die out sufficiently to approximate it as a finite length response.

Here's a plot of the filtered outputs. If you check, you'll see that the outputs for phase_factors 0 and 1 are effectively identical, and the output with phase_factor = 0.5 is very close to the others over the input range 0 < n < 1000, but it has oscillations in the tails, again probably because -NVimp:NVimp isn't long enough for that case.

figure
plot(ny2,vertcat(y2{:}))
xlabel('n');
ylabel('Filtered Output')

function [y2,ny2] = analyze_this(Vfilter,N,x,phase_factor)

Vimp  = impz(Vfilter);
NVimp = numel(Vimp);
figure
nexttile
plot(0:numel(Vimp)-1,Vimp);
xlabel('n')
ylabel('Vimp')

[Vfreq,wfreq] = freqz(Vfilter,N*4,'whole');
Vfreq(end+1)  = Vfreq(1); wfreq(end+1) = 2*pi;
sqrtVfreq     = sqrt(abs(Vfreq)).*exp(1j*angle(Vfreq)*phase_factor);
sqrtVfreqz    = @(w) interp1(wfreq,sqrtVfreq,w);

sqrtVimp = nan(2*NVimp+1,1);
for n = -NVimp:NVimp
    
    ii = n + NVimp + 1;
    sqrtVimp(ii) = trapz(wfreq,sqrtVfreqz(wfreq).*exp(1j*n*wfreq))/2/pi;
    %sqrtVimp(ii) = trapz(wfreq,func(Vfilter,wfreq).*exp(1j*n*wfreq))/2/pi;
    
end
nexttile
plot(-NVimp:NVimp,real(sqrtVimp)),grid
xlabel('n')
ylabel('real(sqrtVImp)')

nexttile
plot(-NVimp:NVimp,imag(sqrtVimp));
xlabel('n')
ylabel('imag(sqrtVimp)');

sqrtVimp = real(sqrtVimp);

x = x(:).';  % ensure row vector
h = sqrtVimp(:).';              % sqrt filter impulse resonse
y2 = conv(fliplr(h),conv(h,x)); % zero-phase, filtered output
Nh = numel(h); % length of impulse response
Nx = numel(x); % length of input

ny2 = (0:numel(y2)-1) - 2*NVimp;  % time of output

[X,w]  = freqz(x,1,N);
Y2     = freqz(y2,1,w).*exp(-1j*ny2(1)*w);

facq = Vfilter.SampleRate;
fvec = w/2/pi*facq;

nexttile
semilogx(1e-3*fvec,db(abs(Y2./X)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
xlabel('Frequency (kHz)');
ylabel('Magnitude')
nexttile
semilogx(1e-3*fvec,180/pi*angle(Y2./X),1e-3*fvec,180/pi*angle(freqz(Vfilter,fvec,facq)));grid
xlabel('Frequency (kHz)');
ylabel('Phase (deg)')

end

[Bruno Luong]

Hmmm so it actually works with phase_factor = 0.5 . I though it did not work when I tried it, possibly I made some misktake.

Thanks again.

[Paul]

You're quite welcome.

Now that the basic structure is working, at least for this example, it seems reasonable that we can further exploit the assumption that the infinite duration impulse response can be approximated as a finite duration response. Dealing with finite duration signals naturally leads to using the DFT/IDFT (implemented by fft/ifft). The code below does exactly this (it only uses ifft, but could use fft) and runs substantially faster. The code also takes advantage of the properties of the DTFT of a non-causal signal that is symmetric around n = 0.

PassbandFrequency   = 1000;
FilterOrder         = 5;
StopbandAttenuation = 60;
facq                = 2e5;
Vfilter = designfilt('lowpassiir', 'FilterOrder', FilterOrder, ...
    'PassbandFrequency', PassbandFrequency, ...
    'StopbandAttenuation', StopbandAttenuation, ...
    'SampleRate', facq);

rng(100);
x = rand(1000,1);

[y3,ny3] = analyze_this2(Vfilter,x);

figure
plot(ny3,y3);
xlabel('n')
ylabel('Zero-phae filter output')

function [y2,ny2] = analyze_this2(Vfilter,x)

% Compute the impulse response of the filter. Assume that Matlab does a
% good job of selecting NVimp such that the Vimp(end) is small enough that
% the remaining elements can be ignored. I have found in other problems
% that sometimes NVimp is too small. Caveat emptor.
Vimp  = impz(Vfilter);
NVimp = numel(Vimp);
figure
nexttile
plot(0:numel(Vimp)-1,Vimp);
xlabel('n')
ylabel('Vimp')

% Compute samples from one period of the frequency response of Vfilter.
% It is assumed that impulse response of the sqrt filter can be adequately
% represented over the range n = -NVimp:Nvimp.
% NVimp can be increased, but symmetry around n=0 is important.
% Use freqz for the computation, though it could just as easily be done with fft. 
N = 2*NVimp + 1;
wfreq = (0:N-1)/N*2*pi;
Vfreq = freqz(Vfilter,wfreq);

% The impulse response of the sqrt filter is:
sqrtVimp = fftshift(ifft(sqrt(abs(Vfreq)),'symmetric'));

nexttile
plot(-NVimp:NVimp,real(sqrtVimp)),grid
xlabel('n')
ylabel('real(sqrtVImp)')

nexttile
plot(-NVimp:NVimp,imag(sqrtVimp));
xlabel('n')
ylabel('imag(sqrtVimp)');

%sqrtVimp = real(sqrtVimp);

x = x(:).';                     % ensure row vector
h = sqrtVimp(:).';              % sqrt filter impulse response

% Maybe it would be faster to use fftconv instead of conv?
y2 = conv(fliplr(h),conv(h,x)); % zero-phase, filtered output
ny2 = (0:numel(y2)-1) - 2*NVimp;  % time of output

% Transform of input and output
[X,w]  = freqz(x,1,N);
Y2     = freqz(y2,1,w).*exp(-1j*ny2(1)*w);

facq = Vfilter.SampleRate;
fvec = w/2/pi*facq;

nexttile
semilogx(1e-3*fvec,db(abs(Y2./X)),1e-3*fvec,db(abs(freqz(Vfilter,fvec,facq))));grid
xlabel('Frequency (kHz)');
ylabel('Magnitude')
nexttile
semilogx(1e-3*fvec,180/pi*angle(Y2./X),1e-3*fvec,180/pi*angle(freqz(Vfilter,fvec,facq)));grid
xlabel('Frequency (kHz)');
ylabel('Phase (deg)')

end

[Bruno Luong]

For better view of non-causal FIR

h2 =  conv(fliplr(h),h);
y2 = conv(h2,x)); % zero-phase, filtered output

[Paul]

In the current incarnation of the code (i.e., equivalent to phase_factor = 0), h is mathematically symmetric, so there's no need for the fliplr, i.e., just do h2 = conv(h,h). Also, there might be slight benefit in accuracy by enforicng that symmetry numerically with something like

% k = (numel(h)-1)/2; h = [h(1:k) h(k+1) fliplr(h(1:k))];

Or maybe both sides around the middle should be averaged.

But, more importantly, that definiton of h2 suggests that we should be able to compute h2 directly. Specifically, let V(w) be the DTFT of the impulse response of Vfilter. Then, I think that it should be true that H2(w) = abs(V(w)). Then, the FIR approximation of h2 should be found from the ifft of frequency domain samples of H2(w).

We currently have (with the above-mentioned modifications)

PassbandFrequency   = 1000;
FilterOrder         = 5;
StopbandAttenuation = 60;
facq                = 2e5;
Vfilter = designfilt('lowpassiir', 'FilterOrder', FilterOrder, ...
    'PassbandFrequency', PassbandFrequency, ...
    'StopbandAttenuation', StopbandAttenuation, ...
    'SampleRate', facq);

Vimp  = impz(Vfilter);
NVimp = numel(Vimp);

N = 2*NVimp + 1;
wfreq = (0:N-1)/N*2*pi;
Vfreq = freqz(Vfilter,wfreq);
sqrtVimp = fftshift(ifft(sqrt(abs(Vfreq)),'symmetric'));
h = sqrtVimp(:).';

k = (numel(h)-1)/2;
h = [h(1:k) h(k+1) fliplr(h(1:k))];
h2 = conv(h,h); % no fliplr
k = (numel(h2)-1)/2;
h2 = [h2(1:k) h2(k+1) fliplr(h2(1:k))];
nh2 = -2*NVimp:2*NVimp;

The proposed approach would be

N = 4*NVimp + 1;
wfreq = (0:N-1)/N*2*pi;
Vfreq = freqz(Vfilter,wfreq);
h2new = fftshift(ifft((abs(Vfreq)),'symmetric'));  % no sqrt
h2new = [h2new(1:k) h2new(k+1) fliplr(h2new(1:k))];

At first glance h2 and h2new look the same

figure
plot(nh2,h2,nh2,h2new)

But they are different enough to give some cause for concern

figure
plot(nh2,h2-h2new)

I don't know if that large of a difference is due to 1) an error in my math, 2) the effect of different assumptions in arriving at h2 vs. h2new, or 3) numerical errors that build up differently between the two methods.

My gut feeling is that (2) is the issue and h2new is the better approach because it starts from the fundamental definition of H2(w) and then makes the FIR approximation, whereas h2 starts from sqrt(H2(w)), makes the FIR approximation, and then exacerbates that approximation by conv(h,h). But that's just my feeling.

If one has a LOT of time on one's hands, I suppose both approaches could be implemented using Symbolic Math Toolbox and see if they really are the same and (3) is the issue. I'm not sure this is a serious suggestion.

I did verify that using h2new does yield very clean mag and phase plots of Y2./X.

[Paul]

@Bruno Luong

Just to tie this off (I hope). I wish I had thought about this problem some more before answering, but having gone through this process the whole situation is much more clear to me.

Let h[n] be a discrete-time signal and H(w) be its DTFT. If h[n] is real and even-symmetric then so will be H(w).

In the FIR example, the Vfilter is a Type-1 filter, which has an even order N, and so its impulse response, v[n], is symmetric around N/2. Let h[n] = v[n+N/2], i.e., h[n] is v[n] shifted to the left by N/2. Then we'll have abs(H(w)) = abs(V(w)) (because time shifts don't affect magnitude of the DTFT) and angle(H(w)) = 0 (due to even symmetry of h[n]). Therefore filtering with h[n] yields the desired zero-phase response, which is basically what we did by using filter with v[n] and shifting the output to the left by N/2 (and zero-padding the input sequence to filter to make sure we collected all of the non-zero filtered output). However, this approach won't work for Type 2-4 FIR filters, because of odd order and/or v[n] not having the desired symmetry. However, for Type 2-4 FIR filters, I'm pretty sure we could use the same approach as for the general IIR case to get a very reasonable result.

For the IIR example, as hopefully is clear from the answer above with h2new, which I'll now call h[n] for simplicity, we're doing the exact same thing. We are defining H(w) = abs(V(w)), and then "finding" the corresponding, real, even-symmetric h[n] for filtering. I say "finding" because we are assuming that h[n] can be adequately represented as a finite duration signal, even though it may actually be infinite duration.