How to concatenate IIR filter (Butterworth) with a Matlab function and evaluate the integral

Question

I have a function of $j\Omega$ and $v$ (velocity), i.e.:

$$A(j\Omega,v)=\frac{(j\Omega)^2}{p + \frac{j\Omega}{v}}$$

where $p$ is the pole.

I want know how $A(j\Omega,v)$ modify its behavior, when I insert another cascade filter $F(j\Omega)$.

In this case I would like study the power of this system via Parseval's theorem:

$$P(v) = \frac{1}{2 \pi} \int_{-\infty}^{+\infty} \big|A(j \Omega,v) F(j\Omega)\big|^2 \text{d} \Omega$$

Until yesterday, I have used always the same $F(j\Omega)$, so the Matlab code was long, but simple:

acc = @(w,v) ...
    abs( ...
    ((1i.*w).^2).* ...
    (1./(.1+1i.*w./v)) .*...
          ((0.0612.*(1i.*w)+1)./ ...
          (3.375e-06.*(1i.*w).^4+0.0001057.*(1i.*w).^3+0.0167.* ...
          (1i.*w).^2+0.0612.*(1i.*w)+1)) ...
    ).^2;


q = zeros(v_min, v_max);

for v = v_min:v_max;
    q(v) = integral(@(w)acc(w,v), 0, Inf)./pi;
end

Where the second part of acc is the filter.

This code works and everything is OK.

Now, I do not want put that filter, but several kind of IIR Butterworth low-pass filters.

1. A series, with $F_c = constant$, but with order N = 1:10.
2. A second series, with $N = constant$, and Fc = [1 2 4 8 16 32 64 128]Hz.

I know that the code for generate these filters is:

h  = fdesign.lowpass('N,F3dB', N, Fc, Fs);
Hd = design(h, 'butter');

but my questions is: how concatenate this filter with a Matlab function, i.e.:

A = @(w,v)(((1i.*w).^2).*(1./(.1+1i.*w./v)));

for evaluate the integral?

Answer 1

The functions fdesign and design are for digital filters. Use butter instead to obtain the filter coefficients of an analog butterworth filter. Then use the function freqs to evaluate the (Laplacian) transfer function of that filter at given frequencies.

To use integral I would first create a function that calculates the transfer function of the concatenation of filters A and B. For example:

function H = filter_cascade(w, v)
    % transfer function of A
    H_A = ((1i.*w).^2)./(.1+1i.*w./v);

    % the other transfer function
    [b,a] = butter(3, 0.3);
    H_B = freqs(b, a, w);

    % The transfer function of the cascade of A and B
    H = H_A .* H_B;

    % Return the magnitude square of the transfer function
    H = H.*conj(H);
end

Then pass a handle to this function to integral like this:

v = 0.1;
q = integral(@(w)filter_cascade(w,v), 0, 1)./pi;

Note that I have changed the upper limit of the integral to one.