Filtri Ideali e Reali

Introduzione

La teoria dei Segnali e un Circuito in regime sinusoidale si incontrano al bar, ma non si salutano.

Tra i sistemi fisici ed in particolare sotto i cosiddetti sistemi LTI (lineari tempo-invariati), è opportuno porre l’accento sui filtri ideali. Un filtro è un sistema LTI selettivo nel dominio delle frequenze.

I filtri ideali sono non causali e sono il riferimento a cui approssimare filtri reali, costituiti da circuiti elettrici (resistori, induttori e condensatori). Dunque, grazie alla teoria dei circuiti, è possibile realizzare e studiare fisicamente i filtri.

Le quattro categorie principali di filtri sono:

1. LP: passa-basso (low-pass)
2. HP: passa-alto (high-pass)
3. BP: passa-banda (band-pass)
4. BR: elimina-banda (band-rejection)

---

La banda è l’intervallo di frequenze in cui il segnale assume valori maggiori di zero. Si considera il solo semiasse positivo.

Decibel

Le scale logaritmiche aiutano nella visualizzazione di grandezze fisiche molto estese su uno (o più) assi. Dalla teoria del circuiti è noto che la potenza elettrica è data dal prodotto tra tensione e corrente: $p=vi$.

Date due potenze $P_1$ e $P_2$, il rapporto di potenze R è:

$$R = \frac{P_2}{P_1} = \bigg(\frac{V_2}{V_1}\bigg)^2 = \bigg(\frac{I_2}{I_1}\bigg)^2$$

Il rapporto di potenze espresso in decibel è:

$$R_{dB} = 10 \log_{10}\bigg(\frac{P_2}{P_1}\bigg) = 10 \log_{10}\bigg(\frac{V_2}{V_1}\bigg)^2 = 20 \log_{10}\bigg(\frac{V_2}{V_1}\bigg)$$

R_dB rappresenta il valore di R espresso in decibel.

---

Il decibel (dB) è la decima parte del bel (simbolo B):

$$10 [dB] = 1 [B]$$

Il decibel viene adottato per evidenziare rapporti di potenza.

Se si vuole ottenere un valore assoluto, si sostituisce P₁ con il riferimento di 1 milli-watt e si ottiene il dB_m:

$$dB_{m} = 10 \log_{10} \bigg(\frac{P_2}{1\;[mW]}\bigg)$$

Il dB_m può anche essere indicato come dB_mW per esplicitare ancora meglio che si tratta di un rapporto di potenze.

A livello teorico, la definizione del dB per le tensioni (calcolate sulla stessa resistenza) è:

$$dB = 20 \log_{10} \bigg(\frac{V_2}{V_1}\bigg)$$

Si ottiene, come nel caso precedente, una grandezza relativa. Se si vuole ottenere una grandezza in valore assoluto, allora si deve porre a denominatore 1 milli-volt:

$$dB_{mV} = 20 \log_{10} \bigg(\frac{V_2}{1\;[mV]}\bigg)$$

---

Data la funzione $X(f)$ in ingresso al filtro e $Y(f)$ in uscita dal filtro. $H(f)$ è la risposta in frequenza del filtro ed è il rapporto tra grandezze fisiche che rappresentano delle ampiezze. La risposta in ampiezza di un filtro, ovvero il modulo di $H(f)$, può essere scritta come:

$$A(f) = |H(f)| = \frac{|Y(f)|}{|X(f)|} \implies A(f)_{dB} = 20 \log_{10} A(f)$$

$R=V_2^2 / V_1^2$	$V_2/V_1$	$R_{dB}$
$10$	$\sqrt{10}$	$10$
$100$	$10$	$20$
$1000$	$10\sqrt{10}$	$30$
$1$	$1$	$0$
$0.1$	$\sqrt{0.1}$	$-10$
$0.01$	$0.1$	$-20$
$2$	$\sqrt{2}$	$3$
$1/2$	$1/\sqrt{2}$	$-3$

Dalla tabella soprastante si nota che la banda a $-3 dB$ del segnale $H(f)$ si ottiene $A(f)$.

Grafici di Filtri Fondamentali

In questo articolo sono presenti simulazione effettuate con Simulink e codice scritto in MATLAB.

%%%
% sample usage:
% frequency_response(10, 2, 1, 4)
%
% parameters:
% #1: horizontal plots limit (both positive and negative)
% #2: band value
% #3: lower bound for band-pass and band-rejection plots
% #4: upper bound for band-pass and band-rejection plots
%%%

function frequency_response(f_lim, B, fL, fH)

    % fprintf('»»» number of arguments: %d «««\n', nargin);

    if nargin ~= 4
        fprintf(['»»» USAGE: ideal_filters(f_lim, B, fL, fH)\n' ...
            '\tf_lim:\thorizontal plots limit (both positive and negative)\n' ...
            '\tB:\tband value\n' ...
            '\tfL:\tlower bound for band-pass and band-rejection plots\n' ...
            '\tfH:\tupper bound for band-pass and band-rejection plots\n' ...
            '«««\n']);
        return;
    end

    % take the absolute value of the first parameter and set as horizontal axes
    f_lim = abs(f_lim);
    f = linspace(-f_lim, f_lim, 1001);

    % define the functions
    H_lp = low_pass(f, B); % rect(f ./(2 * B));
    H_hp = high_pass(f, B);
    H_bp = band_pass(f, fL, fH);
    H_br = band_rejection(f, fL, fH);

    % ---

    figure('Name', 'Ideal Filters');
    plots = [H_lp; H_hp; H_bp; H_br];
    titles = ["Low Pass" "High Pass" "Band Pass" "Band Rejection"];

    % the use of i1 as loop indexes could cause some issuces since
    % it is mainly used to represent the immaginary unit
    for a = 1:4
        % fprintf('»»» i1 = %d', i1)
        subplot(2,2,a)
        plot(f, plots(a,:),'Color','blue','LineWidth',1.5)
        hold on
        plot(f, zeros(size(f)),zeros(size(f)),f,'Color','#808080','LineStyle','--','LineWidth',1.5)
        title(titles(a))
        hold off
        grid on
        ylim([-0.1, 1.1])
    end

end

% low_pass_plot.DataTipTemplate.DataTipRows(1).Label = "-B"
% datatip(lp_plot,-B,0);
% datatip(lp_plot,B,0);

function y = low_pass(x, B)
    y = rect(x ./ (2 * B));
    % den = 2 * B; % denominator
    % y = zeros(1, length(x));
    % y(abs(x./den) < 0.5) = 1;
    % y(abs(x./den) == 0.5) = 0.5;
end

function y = high_pass(f, B)
    y = 1 - low_pass(f, B);
    % den = 2 * B; % denominator
    % y = zeros(1, length(x));
    % y(abs(x./den) > 0.5) = 1;
    % y(abs(x./den) == 0.5) = 0.5;
end

function y = band_pass(f, fL, fH)
    B = fH - fL;
    f0 = (fH + fL) / 2;

    rect1 = rect((f-f0) ./ B);
    rect2 = rect((f+f0) ./ B);
    y = rect1 + rect2;
end

function y = band_rejection(f, fL, fH)
    y = 1 - band_pass(f, fL, fH);
end

function y = rect(x)
    y = zeros(1, length(x));
    y(abs(x) < 0.5) = 1;
    y(abs(x) == 0.5) = 0.5;
end

Le risposte in frequenza dei quattro filtri ideali più noti sono:

%%%
% sample usage:
% impulse_response(10, 2, 1, 4)
%
% parameters:
% #1: horizontal plots limit (both positive and negative)
% #2: band value
% #3: lower bound for band-pass and band-rejection plots
% #4: upper bound for band-pass and band-rejection plots
%%%

function impulse_response(t_lim, B, fL, fH)
    if nargin ~= 4
        fprintf(['»»» USAGE: impulse_response(t_lim, B, fL, fH)\n' ...
            '\tt_lim:\thorizontal plots limit (both positive and negative)\n' ...
            '\tB:\tband value\n' ...
            '\tfL:\tlower bound for band-pass and band-rejection plots\n' ...
            '\tfH:\tupper bound for band-pass and band-rejection plots\n' ...
            '«««\n']);
        return;
    end

    % take the absolute value of the first parameter and set as horizontal axes
    t_lim = abs(t_lim);
    t = linspace(-t_lim, t_lim, 10001);
    B2 = (2 * B)+0.2;

    % define the functions
    h_lp = impulsive_low_pass(t, B);
    h_hp = impulsive_high_pass(t, B);
    h_bp = impulsive_band_pass(t, fL, fH);
    h_br = impulsive_band_rejection(t, fL, fH);

    % ---

    figure('Name', 'Impulsive Responses');
    plots = [h_lp; h_hp; h_bp; h_br];
    titles = ["Low Pass" "High Pass" "Band Pass" "Band Rejection"];

    % create separate plots for each impulse response
    % the outer loop determinates which figure to plot (all togheter or not)
    % the inner loop draws the subplots
    % the use of i1 and j1 as loop indexes could cause some issuces since
    % they are mainly used to represent the immaginary unit
    for a = 1:2
        % fprintf('»»» \n i1 = %d', i1)

        for b = 1:4
            % fprintf('»»» \n j1 = %d', j1)

            if a == 1
                subplot(2,2,b)
            elseif a == 2
                figure('Name', titles(b))
            end

            plot(t, plots(b,:), 'Color', 'blue','LineWidth',1.5)
            hold on
            plot(t, zeros(size(t)), zeros(size(t)), t, 'Color', '#808080', 'LineStyle','--','LineWidth',1.5)
            title(titles(b))
            hold off
            grid on
            ylim([-B2, B2])
        end
    end
end

function y = impulsive_low_pass(x, B)
    y = 2 * B * sinc(2*B .* x);

    % true_h_lp = zeros(1,length(t));
    % two_b = B * 2;
    % true_h_lp(t==0)=1;
    % true_h_lp(t~=0)= two_b * sin(two_b * pi*t(t~=0))./(two_b * pi*t(t~=0));
end

function y = impulsive_high_pass(x, B)
    y = dirac(x) - 2 * B * sinc(2*B .* x);
end

function y = impulsive_band_pass(x, fL, fH)
    B = fH - fL;
    f0 = (fH + fL) / 2;

    y = 2 * B * sinc(B .* x) .* cos(2 * pi * f0 .* x);
end

function y = impulsive_band_rejection(x, fL, fH)
    y = dirac(x) - impulsive_band_pass(x, fL, fH);
end

function y = sinc(x)
    y = zeros(1, length(x));
    y(x==0) = 1;
    y(x~=0) = sin(pi*x(x~=0))./(pi*x(x~=0));
end

Le risposte impulsive (nel dominio temporale) dei quattro filtri ideali più noti sono:

Passa-Basso Ideale

Simbologia: due sinusoidi con quella più in alto barrata. Le frequenze più alte sono tagliate. La sigla LP in pedice significa low pass.

Il parametro B è la banda del filtro passa-basso.

$$H_{\rm LP}(f) = \text{rect}\bigg(\frac{f}{2B}\bigg) \quad ; \quad h_{\rm LP}(t) = 2B \cdot \text{sinc}(2Bt)$$

La banda passante è compresa tra $0$ e B. La banda oscura si estende da B a $\infty$.

Passa-Basso Reale RC

Un filtro passa-basso reale è un circuito formato da un resistore con valore di resistenza R in serie ad un condensatore con capacità elettrica C. Lo schema elettrico potrebbe trarre in inganno facendo sembrare i componenti in parallelo. Il segnale in ingresso $x(t)$ è la tensione ai capi dei due componenti in serie, mentre il segnale di uscita $y(t)$ è la tensione ai capi del solo condensatore. La corrente $i(t)$ scorre nel circuito (idealmente dall’ingresso all’uscita).

Il 1° metodo di risoluzione consiste nel considerare il sistema LTI:

$$y(t) = \mathcal{T} \Big[ x(t)\Big]$$

Al circuito è associato il seguente sistema di equazioni:

$$\sis{ i(t)=C \frac{d}{dt}y(t)\\ \\ x(t)=y(t)+R i(t) }$$

La risposta impulsiva del sistema passa-basso reale è la risposta all’impulso Delta di Dirac. Inoltre, si ricorda che la Delta di Dirac è la derivata della funzione gradino.

$$\Large \xrightarrow{u(t)} \boxed{\frac{d}{dt}} \xrightarrow{\delta(t)} \boxed{\mathcal{T}} \xrightarrow{h(t)}$$

Un sistema LTI a cascata è invertibile, dunque:

$$\Large \xrightarrow{u(t)} \boxed{\mathcal{T}} \xrightarrow{z(t)} \boxed{\frac{d}{dt}} \xrightarrow{h(t)}$$

con $z(t)\neq \delta (t)$.

Allora $z(t)=\mathcal{T}[u(t)]$ rappresenta il transitorio di carica del condensatore:

$$u(t)=\sis{1 \quad t \geq 0 \\ 0\quad t < 0} \implies z(t)=\sis{ 1-e^{t/RC}\quad &t \geq 0 \\ 0 &t < 0 }$$

La risposta impulsiva è la derivata del transitorio di carica:

$$\Large h(t)=\frac{d}{dt}z(t) = \sis{ \frac{1}{RC}\cdot e^{-t/RC}\quad &t \geq 0 \\ 0 &t < 0 }$$

Si nota la notevole somiglianza con la definizione di esponenziale monolatero:

$$\text{monoExp}(t) = e^{-t/T} \cdot u(t)= \sis{e^{-t/T} \quad & t > 0\\ 0 & t \leq 0 }$$

La risposta in frequenza è la trasformata di Fourier della risposta impulsiva:

$$\eq{ H(f) &= \int_{-\infty}^{+\infty}x(t)e^{-j2\pi ft}dt =\\ &= \int_{0}^{+\infty} \frac{1}{RC} e^{-t/(RC)} e^{-j2\pi ft}dt =\\ &= \frac{1}{RC}\int_{0}^{+\infty} e^{-t(1/(RC)+j2\pi f)}dt =\\ &= \frac{1}{RC}\lim_{c \to +\infty}\int_{0}^c e^{-t(1/(RC)+j2\pi f)}dt=\\ &= \frac{1}{RC}\frac{RC}{1+j2\pi fRC} =\\ &= \frac{1}{1+j2\pi fRC} }$$

La risposta in ampiezza è il modulo della risposta in frequenza:

$$A(f) = \frac{1}{\sqrt{1+(2\pi f RC)^2}}$$

La risposta in fase è l’argomento della risposta in frequenza:

$$\Theta(f) = \arg \bigg(H(f)\bigg) = -\text{arctan}(2\pi f RC)$$

---

Il 2° metodo di risoluzione consiste nell’utilizzo del metodo dei fasori, il quale prevede il passaggio ai numeri complessi per rappresentare circuiti elettrici che operano in regime sinusoidale. L’analisi circuitale in regime sinusoidale prevede che ad ogni elemento sia associato un fasore, un valore di impedenza e la corrispondente ammettenza.

L’articolo Analisi in Regime Sinusoidale tratta in modo molto approfondito questo argomento. Si può usare la tabella che segue per ottenere l’impedenza relativa ad un componente:

Elemento	Tempo	Fasore	Impedenza	Ammettenza
R: Resistore	$v=Ri$	$V=RI$	$Z = R$	$Y=1/R$
L: Induttore	$v = L \frac{di}{dt}$	$V = j \omega LI$	$Z = j \omega L$	$Y=1/(j \omega L)$
C: Condensatore	$i = C \frac{dv}{dt}$	$V = I/(j \omega C)$	$Z = 1/(j \omega C)$	$Y = j \omega C$

Da notare che la pulsazione è $\omega = 2\pi f$. Nelle formule che seguono, si esplicita la $f$ in quanto variabile indipendente dei grafici $H(f)$, $A(f)$ e $\Theta(f)$.

Si applica un partitore tra i valori di impedenza della resistenza e del condensatore. Dato il fasore V_X della funzione in ingresso $x(t)$. Il fasore V_Y della tensione in uscita $y(t)$ è dato da:

$$V_Y = V_X \cdot \frac{Z_C}{Z_R + Z_C}$$

La risposta in frequenza è il rapporto tra il fasore della tensione in uscita e quello della tensione in ingresso:

$$H(f) = \frac{V_Y}{V_X} = \frac{Z_C}{Z_R + Z_C} = \frac{1}{1+j2\pi fRC}$$

La risposta in ampiezza è il modulo della risposta in frequenza:

$$A(f) = \frac{1}{\sqrt{1+(2\pi f RC)^2}}$$

La risposta in fase è l’argomento della risposta in frequenza:

$$\Theta(f) = \arg \bigg(H(f)\bigg) = -\text{arctan}(2\pi f RC)$$

%%%
% sample usage:
% real_low_pass_RC(1, 0.1)
%
% parameters:
% #1: resistor value in Ohm
% #2: capacitor value in Farad
%%%

function real_low_pass_RC(R, C)
  f = linspace(-40, 40, 10001);
    f_size = zeros(size(f));
    w = 2*pi.*f;
    Z_R = R;
    Z_C = 1./(1j*w*C);
    Hf = Z_C ./ (Z_R + Z_C);

    % amplitude response: A(f) = |H(f)|
    Af = 1 ./ sqrt(1 + (2 * pi * R * C .* f) .^ 2);

    % the filter is a low pass: f0 = 0;
    Af0 = 1 / sqrt(2);

    % B_(-3) --> band value at -3 dB
    B3 = 1 / (2 * pi * R * C);

    figure(Name="Real Low Pass");
    tiledlayout(1, 2)
    nexttile
    plot(f, Af, 'LineWidth', 2.5, 'Color', 'blue')
    hold on
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    hold on
    yline(Af0, 'Color', 'red', 'LineStyle','-.')
    hold on
    xline(B3, 'Color', 'red', 'LineStyle','-.')
    hold off
    grid on
    ylim([-0.1, 1.1])
    xlim([-20 20])
    ax = gca;
    chart = ax.Children(5);
    datatip(chart,B3,Af0);
    chart = ax.Children(4);
    datatip(chart,B3,0);
    chart = ax.Children(3);
    datatip(chart,0,Af0,"Location","southwest");
    title('Real LP - Amplitude Response', FontSize=20)

    % ---

    % risposta in fase
    nexttile
    Fasef = angle(Hf);
    plot(f, Fasef, 'LineWidth', 2.5, 'Color', 'red')
    hold on
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    hold off
    grid on
    ylim([-2 2])
    xlim([-20 20])
    title('Real LP - Phase Response', FontSize=20)
end

Passa-Basso Reale LR

Un filtro passa-basso reale, a differenza di quanto visto sopra, può anche essere formato ad un induttore con valore di induttanza L in serie ad un resistore con valore di resistenza R. Il segnale in ingresso $x(t)$ è la tensione ai capi dei due componenti in serie, mentre il segnale di uscita $y(t)$ è la tensione ai capi del solo resistore.

Utilizzando il metodo dei fasori, si ottiene il fasore della tensione in uscita:

$$V_Y = V_X \cdot \frac{Z_R}{Z_R + Z_L}$$

Nel dominio delle frequenze si può scrivere:

$$Y(f) = X(f) \cdot \frac{Z_R}{Z_R + Z_L} = X(f) \cdot \frac{R}{R + j 2 \pi fL}$$

La risposta in frequenza è il rapporto tra il fasore della tensione in uscita e quello della tensione in ingresso:

$$\eq{ H(f)&=\frac{V_Y}{V_X}=\frac{Y(f)}{X(f)}=\frac{Z_R}{Z_R + Z_L} =\\ &=\frac{R}{R+j2\pi fL}=\frac{1}{1+(j2\pi fL)/R}=\\ &=\frac{1}{1+j2\pi f \frac{L}{R}} }$$

La risposta in ampiezza è il modulo della risposta in frequenza:

$$A(f) = |H(f)| = \frac{R^2}{\sqrt{R^2 + (2\pi fL)^2}} = \frac{1}{\sqrt{1+\Big(2\pi f \frac{L}{R}\Big)^2}}$$

Studiando la risposta in ampiezza al limite, si ricava che:

$$|H(f)| \xrightarrow{f \to 0} 1 \quad ; \quad |H(f)| \xrightarrow{f \to \infty} 0\\$$

La risposta in ampiezza trova il proprio massimo per $f_0 = 0$, quindi si tratta di un filtro passa-basso.

$f_{-3}$ si ottiene grazie alla relazione che segue:

$$\frac{|H(f_{-3})|}{|H(f_0)|} = \frac{1}{\sqrt{2}}$$ $$\eq{ |H(f_{-3})| &= \frac{1}{\sqrt{2}} |H(f_0)| \\ \implies & \frac{1}{\sqrt{1+\Big(2\pi f_{-3} \frac{L}{R}\Big)^2}}=\frac{1}{\sqrt{2}}\\ \implies & 1+\Big(2\pi f_{-3} \frac{L}{R}\Big)^2 = 2\\ \implies & f_{-3}^2\Big(2\pi\frac{L}{R}\Big)^2 = 1 \\ \implies & f_{-3}^2 = \frac{R^2}{4\pi^2L^2} \\ \implies & f_{-3}= \pm \frac{R}{2\pi L} \\ }$$

La banda a $-3dB$ del filtro è il valore positivo (tra i due possibili) di $f_{-3}$:

$$B_{-3 dB} \equiv f_{-3}= \frac{R}{2\pi L}$$

La risposta in fase è l’argomento della risposta in frequenza:

$$\Theta(f) = \arg \bigg(H(f)\bigg) = -\text{arctan}\bigg(\frac{2\pi f L}{R}\bigg)$$

%%%
% sample usage:
% real_low_pass_LR(1, 12.57)
%
% parameters:
% #1: inductor value in Henry
% #2: resistor value in Ohm
%%%

function real_low_pass_LR(L, R)
    f = linspace(-40, 40, 10001);
    f_size = zeros(size(f));
    % pulsazione (omega)
    w = 2*pi.*f;
    Z_R = R;
    Z_L = 1j*w*L;
    Hf = Z_R ./ (Z_R + Z_L);
    % amplitude response: A(f) = |H(f)|
    Af = 1 ./ sqrt(1 + (2 * pi * (L/R) .* f) .^ 2);
    Af0 = 1 / sqrt(2);

    % B_(-3) --> band value at -3 dB
    B3 = R / (2 * pi * L);
    fprintf("»»» valore di banda B = %f\n", B3);

    figure(Name="Real Low Pass - LR");
    tiledlayout(1, 2)
    nexttile
    plot(f, Af, 'LineWidth', 2.5, 'Color', 'blue')
    % ...
    % risposta in fase
    nexttile
    Fasef = angle(Hf);
    plot(f, Fasef, 'LineWidth', 2.5, 'Color', 'red')
    % ...
end

Calcolo dell’uscita dato l’ingresso

$$\xrightarrow{x(t)} \boxed{\rm Low \; Pass} \xrightarrow{y(t)}$$

Dato un ingresso nella forma:

$$x(t) = \alpha \cos(2\pi f_0 t)$$

N.B. La $f_0$ a cui si fa riferimento in questo paragrafo non è quella del paragrafo precedente.

Allora l’uscita è data da:

$$y(t) = \alpha \Big|H(f_0)\Big| \cos\Big(2\pi f_0 t+\arg\big(H(f_0)\big)\Big)$$

Esempio:

Sia $x(t) = 8 \cos(2\pi f_0 t)$ con $2\pi f_0=10^3 s^{-1}$ e $L/R = \sqrt{3}\cdot 10^{-3}s$.

Allora $\alpha = 8$:

$$2\pi f_0=10^3 \implies f_0 = \frac{10^3}{2\pi}\\ \frac{L}{R} = \frac{\sqrt{3}}{10^3}$$

La risposta in ampiezza è:

$$\eq{ \Big|H(f_0)\Big| &= \bigg|H\bigg(\frac{10^3}{2\pi}\bigg)\bigg|=\\ &= \frac{1}{\sqrt{1+\Big(2\pi \frac{10^3}{2\pi} \frac{\sqrt{3}}{10^3}\Big)^2}}=\\ &= \frac{1}{\sqrt{1+(\sqrt{3})^2}} = \frac{1}{\sqrt{4}}=\frac{1}{2} }$$

La risposta in fase è:

$$\eq{ \Theta(f_0) &= -\text{arctan}\bigg(\frac{2\pi f_0 L}{R}\bigg) =\\ &= -\text{arctan}\bigg(2\pi \frac{10^3}{2\pi} \frac{\sqrt{3}}{10^3}\bigg) =\\ &= -\text{arctan}\Big(\sqrt{3}\Big) =\\ &= -\frac{\pi}{3}\\ }$$

Allora l’uscita è data da:

$$\eq{ y(t) &= \alpha \Big|H(f_0)\Big| \cos\Big(2\pi f_0 t+\arg\big(H(f_0)\big)\Big)=\\ &= 8 \frac{1}{2}\cos\bigg(2\pi \frac{10^3}{2\pi} t-\frac{\pi}{3}\bigg)=\\ &= 4\cos\bigg(10^3 t-\frac{\pi}{3}\bigg) }$$

Passa-Alto Ideale

Simbologia: due sinusoidi con quella più in basso barrata. Le frequenze più basse sono tagliate. La sigla HP in pedice significa high pass.

$$H_{\rm HP}(f) = 1-rect\bigg(\frac{f}{2B}\bigg) \quad ; \quad h_{\rm HP}(t) = \delta(t) -2B \cdot sinc(2Bt)$$

La banda passante si estende all’infinito. La banda oscura è compresa tra $0$ e B.

Passa-Alto Reale CR

Un filtro passa-alto reale è un circuito formato da un condensatore con capacità elettrica C in serie ad un resistore con valore di resistenza R. Lo schema elettrico potrebbe trarre in inganno facendo sembrare i componenti in parallelo. Il segnale in ingresso $x(t)$ è la tensione ai capi dei due componenti in serie, mentre il segnale di uscita $y(t)$ è la tensione ai capi del solo resistore. La corrente $i(t)$ scorre nel circuito (idealmente dall’ingresso all’uscita).

Utilizzando il metodo dei fasori come spiegato per il passa basso, si ottiene il fasore della tensione in uscita. In questo caso, nella formula del partitore, l’impedenza del resistore e l’impedenza del condensatore sono invertite.

$$V_Y = V_X \cdot \frac{Z_R}{Z_R + Z_C}$$

La risposta in frequenza è il rapporto tra il fasore della tensione in uscita e quello della tensione in ingresso:

$$H(f) = \frac{V_Y}{V_X} = \frac{Z_R}{Z_R + Z_C} = \frac{1}{1-j2\pi fRC}$$

La risposta in ampiezza è il modulo della risposta in frequenza:

$$\Large A(f) = |H(f)| = \frac{1}{\sqrt{1+\frac{1}{(2\pi f RC)^2}}}$$

Studiando la risposta in ampiezza per $f \to 0$ e poi per $f \to \infty$ si ottengono i valori:

$$\eq{ \widehat{f}_0 &=0\implies |X(\widehat{f}_0)| = 0\\ \widetilde{f}_0 &=\infty\implies |X(\widetilde{f}_0)| = 1\\ }$$

Allora si ottiene:

$$\max|X(f)| = |X(f_0)| = 1\implies f_0 = \widetilde{f}_0 = \infty$$

La risposta in fase è l’argomento della risposta in frequenza:

$$\Theta(f) = \arg \bigg(H(f)\bigg) = \text{arctan}\bigg(\frac{1}{2\pi fRC}\bigg)$$

%%%
% sample usage:
% real_high_pass(1, 0.1)
%
% parameters:
% #1: resistor value in Ohm
% #2: capacitor value in Farad
%%%

function real_high_pass(R, C)
    if nargin ~= 2
        fprintf(['»»» USAGE: real_high_pass(R, C)\n' ...
            '\tR:\tresistor value [Ohm]\n' ...
            '\tC:\tcapacitor value [Farad]\n' ...
            '«««\n']);
        return;
    end

    f = linspace(-40, 40, 10001);
    f_size = zeros(size(f));
    w = 2*pi.*f; % pulsazione (omega)
    Z_R = R;
    Z_C = 1./(1j*w*C);
    Hf = Z_R ./ (Z_R + Z_C);

    % amplitude response: A(f) = |H(f)|
    Af = 1 ./ sqrt(1 +(1./(2 * pi * R * C .* f) .^ 2));

    Af0 = 1 / sqrt(2);

    % B_(-3) --> band value at -3 dB
    B3 = 1 / (2 * pi * R * C);
    fprintf("»»» valore di banda B = %f\n", B3);
    figure(Name="Real High Pass");
    tiledlayout(1, 2)
    nexttile
    plot(f, Af, 'LineWidth', 2.5, 'Color', 'blue')
    hold on
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    hold on
    yline(Af0, 'Color', 'red', 'LineStyle','-.')
    hold on
    xline(B3, 'Color', 'red', 'LineStyle','-.')
    hold off
    grid on
    ylim([-0.1, 1.1])
    xlim([-20 20])
    ax = gca;
    chart = ax.Children(5);
    datatip(chart,B3,Af0);
    chart = ax.Children(4);
    datatip(chart,B3,0);
    chart = ax.Children(3);
    datatip(chart,0,Af0,"Location","southwest");
    title('Real HP - Amplitude Response', FontSize=20)

    % ---

    % risposta in fase
    nexttile
    Fasef = angle(Hf);
    plot(f, Fasef, 'LineWidth', 2.5, 'Color', 'red')
    hold on
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    hold off
    grid on

    ylim([-2 2])
    xlim([-20 20])
    title('Real HP - Phase Response', FontSize=20)
end

Passa-Alto RL

Un filtro passa-alto reale di tipo RL è un circuito formato da un resistore con valore di resistenza R in serie ad un induttore con valore di induttanza L. Lo schema elettrico potrebbe trarre in inganno facendo sembrare i componenti in parallelo. Il segnale in ingresso $x(t)$ è la tensione ai capi dei due componenti in serie, mentre il segnale di uscita $y(t)$ è la tensione ai capi del solo induttore. La corrente $i(t)$ scorre nel circuito (idealmente dall’ingresso all’uscita).

Utilizzando il metodo dei fasori come spiegato per il passa basso, si ottiene il fasore della tensione in uscita.

$$V_Y = V_X \cdot \frac{Z_L}{Z_R + Z_L}$$

La risposta in frequenza è il rapporto tra il fasore della tensione in uscita e quello della tensione in ingresso:

$$H(f) = \frac{V_Y}{V_X} = \frac{Z_L}{Z_R + Z_L} = \frac{j2\pi fL}{R+j2\pi fL} = \frac{1}{1+R/(j2\pi fL)}$$

La risposta in ampiezza è il modulo della risposta in frequenza:

$$\Large A(f) = \frac{1}{\sqrt{1+\Big(\frac{R}{2\pi fL}\Big)^2}}$$

La risposta in fase è l’argomento della risposta in frequenza:

$$\Theta(f) = \arg \bigg(H(f)\bigg) = \text{arctan}\bigg(\frac{R}{2\pi fL}\bigg)$$

%%%
% sample usage:
% real_high_pass_RL(50, 1)
%
% parameters:
% #1: resistor value in Ohm
% #2: inductor value in Henry
%%%

function real_high_pass_RL(R, L)
    if nargin ~= 2
        fprintf(['»»» USAGE: real_high_pass_RL(R, L)\n' ...
            '\tR:\tresistor value [Ohm]\n' ...
            '\tL:\tinductor value [Henry]\n' ...
            '«««\n']);
        return;
    end

    f = linspace(-40, 40, 10001);
    f_size = zeros(size(f));
    w = 2*pi.*f; % pulsazione (omega)
    Z_R = R;
    Z_L = 1j*w*L;
    Hf = Z_L ./ (Z_R + Z_L);

    % amplitude response: A(f) = |H(f)|
    Af = 1./(sqrt(1+(R./(2*pi*L .*f)).^2));

    Af0 = 1 / sqrt(2);

    % B_(-3) --> band value at -3 dB
    B3 = R / (2 * pi * L);
    fprintf("»»» valore di banda B = %f\n", B3);
    figure(Name="Real High Pass");
    tiledlayout(1, 2)
    nexttile
    plot(f, Af, 'LineWidth', 2.5, 'Color', 'blue')
    hold on
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    hold on
    yline(Af0, 'Color', 'red', 'LineStyle','-.')
    hold on
    xline(B3, 'Color', 'red', 'LineStyle','-.')
    hold off
    grid on
    ylim([-0.1, 1.1])
    xlim([-20 20])
    ax = gca;
    chart = ax.Children(5);
    datatip(chart,B3,Af0);
    chart = ax.Children(4);
    datatip(chart,B3,0);
    chart = ax.Children(3);
    datatip(chart,0,Af0,"Location","southwest");
    title('Real HP - Amplitude Response', FontSize=20)

    % ---

    % risposta in fase
    nexttile
    Fasef = angle(Hf);
    plot(f, Fasef, 'LineWidth', 2.5, 'Color', 'red')
    hold on
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    hold off
    grid on

    ylim([-2 2])
    xlim([-20 20])
    title('Real HP - Phase Response', FontSize=20)
end

Passa-banda Ideale

Simbologia: tre sinusoidi con quelle esterne barrate. Le frequenze tagliate sono quelle più basse di $f_L$ (low) e più alte di $f_H$ (high). La sigla BP in pedice significa band pass.

La banda passante è compresa tra $f_L$ e $f_H$, si può quindi affermare che $B = f_H-f_L$.

La frequenza di centro-banda è:

$$f_0 = \frac{f_L + f_H}{2}$$

Il fattore di selettività (o qualità) è: $Q = f_0/B$

$$\eq{ H_{\rm BP} &= rect\bigg(\frac{f-f_0}{B}\bigg) + rect\bigg(\frac{f+f_0}{B}\bigg) \\ h_{\rm BP} &= 2B \cdot sinc(Bt) \cdot \cos(2\pi f_0 t) }$$

Passa-Banda Reale RLC

Un bipolo è una rete elettrica caratterizzata da due soli terminali (o poli) di accesso. Può essere costituito da uno o più elementi circuitali connessi fra loro.

Un filtro passa-banda reale è un circuito formato da un resistore con valore di resistenza R in serie ad un bipolo costituito da un condensatore con capacità elettrica C in parallelo ad un induttore con induttanza L. Il segnale in ingresso $x(t)$ è la tensione ai capi dei due componenti in serie (resistore e bipolo), mentre il segnale di uscita $y(t)$ è la tensione ai capi del solo bipolo (condensatore in parallelo ad induttore).

Il metodo dei fasori prevede il passaggio a numeri complessi per rappresentare circuiti elettrici che operano in regime sinusoidale. Si ottengono i fasori:

$$Z_R = R \quad ; \quad Z_L = j 2\pi f L \quad ; \quad Z_C = \frac{1}{j 2\pi f C}$$

Si calcola Z_P, ovvero l’impedenza in parallelo tra Z_L e Z_C:

$$\eq{ Z_P &= Z_L \parallel Z_C = \frac{Z_L \cdot Z_C}{Z_L + Z_C} = \\ &= \frac{j 2\pi f L}{-(2\pi f)^2 LC + 1} }$$

Dato il fasore V_X della funzione in ingresso $x(t)$. Il fasore V_Y della tensione in uscita $y(t)$ è dato da:

$$V_Y = V_X \cdot \frac{Z_P}{Z_R + Z_P}$$

La risposta in frequenza è il rapporto tra il fasore della tensione in uscita e quello della tensione in ingresso:

$$H(f) = \frac{V_Y}{V_X} = \frac{Z_P}{Z_R + Z_P} = \frac{j 2\pi f L}{R[1-(2\pi f)^2 LC ] + j 2\pi f L }$$

La risposta in ampiezza è il modulo della risposta in frequenza:

$$\Large A(f) = |H(f)| = \frac{1}{\sqrt{1+\frac{[1-(2\pi f)^2 LC]^2}{(2\pi f)^2 L^2}\cdot R^2}}$$

La frequenza massima della risposta in ampiezza è:

$$\eq{ f_{\rm max}&=\max[A(f)] \implies 1-(2\pi f)^2 LC = 0 \\ \implies f_{\rm max} &= \pm \frac{1}{2\pi \sqrt{LC}} }$$

La frequenza di centro banda coincide con la frequenza massima, nel semi-piano positivo:

$$f_0 = \frac{1}{2\pi \sqrt{LC}}$$

La frequenza massima (sull’asse orizzontale) interseca il grafico della risposta in ampiezza per $A(f_{\rm max})=1$. Le frequenze $f_L$ ed $f_H$ intersecano il grafico per $A(f)=1/\sqrt{2}$:

$$\frac{A(f_{-3})}{A(f_0)} = \frac{1}{\sqrt{2}}$$

Bisogna fare qualche utile precisazione:

$$\eq{ f_0 \equiv f_{\rm max} \implies & A(f_0)=A(f_{\rm max})=1\\ \implies & \frac{A(f_{-3})}{1} = \frac{1}{\sqrt{2}}\\ \implies & A(f_{-3}) = \frac{1}{\sqrt{2}}\\ }$$

La banda $B_{-3}=-3[dB]$ è data da:

$$\eq{ & f_{-3L} \equiv f_L[dB] \;;\; f_{-3H} \equiv f_H[dB]\\ & B_{-3} = f_{-3H}-f_{-3L}\\ }$$

%%%
% sample usage:
% real_band_pass(5, 0.4, 0.01)
%
% parameters:
% #1: resistor value in Ohm
% #2: inductor value in Henry
% #3: capacitor value in Farad
%%%

function real_band_pass(R, L, C)
    f = linspace(-40, 40, 10001);
    f_size = zeros(size(f));
    w = 2*pi.*f;
    % calcolo di H(f)
    Z_R = R;
    Z_L = 1j*w*L;
    Z_C = 1./(1j*w*C);
    Z_P = (Z_L .* Z_C)./(Z_L + Z_C);

    Hf = Z_P ./ (Z_R + Z_P);

    % amplitude response: A(f) = |H(f)|
    Af = abs(Hf); % amplitude response: A(f) = |H(f)|
    % 1 ./ sqrt(1 + (2 * pi * R * C .* f) .^ 2);

    % center frequency: frequenza di centro banda
    f_0 = 1 / (2*pi*sqrt(L*C)); % x-axis
    % quality factor: Il fattore di selettività (o qualità)
    Q = 1/R * sqrt(L/C);
    % band value
    B = f_0 ./ Q;
    f_low = f_0 - (B/2);
    f_high = f_0 + (B/2);

    fprintf("»»» valore di banda B = %f\n", B);
    fprintf("»»» frequenza di centro banda f0 = %f\n", f_0);
    fprintf("»»» quality factor Q = %f\n", Q);
    fprintf("»»» f_low = %f\n", f_low);
    fprintf("»»» f_high = %f\n", f_high);
    figure('Name', 'Real Low Pass - Amplitude Response');
    plot(f, Af, 'LineWidth', 2.5, 'Color', 'blue')
    hold on
    % draw horizontal and vertical axes
    plot(f, f_size, f_size, f, 'Color', '#808080', 'LineStyle','--')
    % ...
    % risposta in fase
    Fasef = angle(Hf);
    figure('Name', 'Real Low Pass - Phase Response');
    plot(f, Fasef, 'LineWidth', 2.5, 'Color', 'red')
    % ...
end

Elimina-Banda

Simbologia: tre sinusoidi con quella interna barrata. Le frequenze tagliate sono quelle più alte di $f_L$ (low) e più basse di $f_H$ (high). La sigla BR in pedice significa band rejection.

Il valore di banda B e di $f_0$ sono uguali al filtro passa-banda.

$$\eq{ H_{\rm BR}&=1-H_{\rm BP}=1-rect\bigg(\frac{f-f_0}{B}\bigg)-rect\bigg(\frac{f+f_0}{B}\bigg) \\ h_{\rm BR}&=\delta(t)-\Big[2B \cdot sinc(Bt) \cdot \cos(2\pi f_0 t)\Big] }$$

La banda passante è compresa tra $0$ e $f_L$ e si estende da $f_H$ a $\infty$.

Calcolo della Banda

I filtri ideali sono definiti utilizzando il valore B di banda passante. Un filtro reale cerca di approssimare il comportamento di un filtro ideale ed è definito tramite circuito composto da resistenze, condensatori e induttori.

I filtri reali sono costruiti per ottenere un certo valore di banda, che diventa quindi il punto di arrivo dello sviluppo.

Passa-Basso RC e Passa-Alto CR

Il valore B di banda passante di un filtro passa-alto/basso reale è legato ai valori di resistenza R e capacità elettrica C.

La costante di tempo del filtro è: $T=RC$.

Grazie alla natura del filtro, la banda coincide con la frequenza di taglio (cutoff):

$$B = f_0 = \frac{1}{2\pi T}=\frac{1}{2\pi RC}$$

Il calcolo della banda è uguale per entrambi i filtri. La differenza tra passa-alto e passa-basso sta nelle frequenze che possono passare: sopra o sotto il valore di banda.

N.B. La frequenza di taglio viene indicata con la c o lo zero in pedice: $f_c \equiv f_0$.

Passa-Basso LR

Il valore B di banda passante di un filtro passa-basso LR (serie) è legato ai valore di induttanza L e resistenza R. La costante di tempo del filtro è:

$$T = \frac{L}{R}$$

La banda coincide con la frequenza di taglio (cutoff):

$$B = f_c = \frac{1}{2\pi T}=\frac{R}{2\pi L}$$

(Esisterà sicuramente un passa-alto del tipo RL. La sua banda si calcolerà sicuramente nel medesimo modo).

Passa-Banda RLC

In un filtro passa-banda RLC, le correlazioni tra i valori di resistenza R, induttanza L e capacità C forniscono i parametri caratteristici del filtro.

La frequenza di centro banda è:

$$f_0 = \frac{1}{2\pi \sqrt{LC}}$$

Il fattore di selettività (o qualità) è:

$$Q = \frac{1}{R} \sqrt{\frac{L}{C}}$$

La banda passante è:

$$B = \frac{f_0}{Q}$$

Le frequenze di cutoff inferiore e superiore sono:

$$f_L = f_0 - \frac{B}{2} \quad ; \quad f_H = f_0 + \frac{B}{2}$$

Il coefficiente di smorzamento è:

$$\zeta = \frac{1}{2Q} = \frac{R}{2} \sqrt{\frac{C}{L}}$$

La frequenza angolare (natural frequency) è:

$$\omega_0 = \frac{1}{\sqrt{LC}}$$

La costante di tempo si ottiene grazie al coefficiente di smorzamento e la frequenza angolare:

$$T = \frac{1}{\omega_0 \sqrt{1-\zeta^2}}$$

Applicazioni pratiche

I filtri studiati nella prima parte di questo articolo trovano importanti applicazioni “reali” nell’editing di immagini e nella modifica di file audio.

Filtri applicati a file audio

Lo script filtro_audio, che verrà illustrato alla fine del paragrafo, permette di applicare i quattro principali tipi di filtro ad una traccia audio fornita dall’utente. Lo script è stato opportunamente modificato, rispetto alla versione illustrata in classe, per permettere all’utente di selezionare i seguenti parametri:

• filename: file audio di ingresso,
• filter_type: tipo di filtro da applicare al file di ingresso,
• B: dimensione della banda in Hz,
• f0: frequenza centrale in Hz.

Si prenda come esempio di file audio un arpeggio, una serie di accordi suonati con una chitarra classica o la canzone Scar Tissue dei Red Hot Chili Peppers.

Passando come filter_type il valore bp o band-pass, lo script filtro_audio filtra le frequenze al di fuori della banda B centrate nella frequenza $f_0$. Si può verificare che, diminuendo progressivamente il valore di B, l’audio in uscita risulterà sempre più selettivo attorno alla nota musicale relativa alla frequenza f0 passata come parametro. Mantenendo invece B costante e facendo variare f0, grazie al grafico delle trasformate dei segnali, è possibile apprezzare lo scostamento della banda passante.

%% Utilizzo:
% filtro_audio("ScarTissue.wav", 'bp', 3.0e3, 3.0e3)
% filename ~ nome del file da leggere (path, URL, ...)
% filter_type ~ una stringa tra le seguenti:
% 'bp','band-pass', 'br','band-rejection', 'lp','low-pass', 'hp', 'high-pass'
% B ~ dimensione di banda [Hz]
% f0 ~ frequenza centrale [Hz]
function filtro_audio(filename, filter_type, B, f0)
% --- Lettura campioni --- %
info = audioinfo(filename);      % informazioni sul file audio
freqC = 44.1e3;                  % [Hz]: frequenza di campionamento
tempoC = 1 / freqC;              % [s]: tempo di campionamento
durata = info.Duration;          % [s]: durata del segnale audio
numCampioni = freqC * durata;    %
%numCampioni=info.TotalSamples;  % numero totale di campioni
% funziona solo se il file è campionato a 44.1e3, altrimenti dà errore

[xstereo,fc] = audioread(filename,[1,numCampioni]);
x = (xstereo(:,1))'; % solo canale sinistro
% la frequenza di campionamento dipende dal file sorgente, es. CD: 44.1 kHz
fprintf('»»» frequenza di campionamento = %d [Hz]\n', fc);

T = 50/B; % durata risposta impulsiva [s]
tempoFiltro=0:tempoC:T;
tempoNorm = tempoFiltro-T/2;
% --- applica un filtro in base al parametro filter_type --- %
switch filter_type
    case {'bp','band-pass'} % --- Filtro passa-banda ideale --- %
        % più opzioni per uno stesso blocco
        figure(Name="Filtro passa banda")
        g=2*B*sinc(B*tempoNorm).*rectpuls(tempoNorm/T).*cos(2*pi*f0*tempoNorm);
    case {'br','band-rejection'} % --- Filtro elimina-banda ideale --- %
        figure(Name="Filtro elimina banda")
        g=2*B*sinc(B*tempoNorm).*rectpuls(tempoNorm/T).*(cos(2*pi*f0*tempoNorm).^2);
    case {'lp','low-pass'}  % --- Filtro passa-basso --- %
        figure(Name="Filtro passa basso")
        g=2*B*sinc(2*B*tempoNorm).*rectpuls(tempoNorm/T);
    case {'hp', 'high-pass'}% --- Filtro passa-alto --- %
        figure(Name="Filtro passa alto")
        g=2*B*sinc(B*tempoNorm).*rectpuls(tempoNorm/T).*cos(2*pi*f0*tempoNorm);
    otherwise
        warning('Unexpected filter type.')
        return;
end

% --- Uscita filtrata --- %
w = conv(g,x) * tempoC; % convoluzione tra g ed x
w = w(length(g):length(w));

tempo = 0:tempoC:durata-tempoC;
tempoW = tempo(1:length(w))+T/2;

durataTransitorio = 0.0250; % [s]
A = 5; % amplificazione: rende più visibile il segnale filtrato (è un valore empirico)

tiledlayout(1, 2); nexttile;
plot((tempo-durataTransitorio)*1e3, x, 'Color', 'cyan', 'LineWidth', 2.5);
hold on;
plot((tempoW-durataTransitorio)*1e3, A*w, 'Color', 'black', 'LineWidth', 1.5);
grid on;
xlabel('Tempo (ms)', FontSize=14);
ylabel('Segnali temporali', FontSize=14);
legend('x(t)', 'A\cdot{w(t)}', FontSize=12);
axis([2000 2200 -1*max(abs(w)) max(abs(w)) ]);
titoloGrafico=sprintf('Filtraggio %s con banda B = %.2f [Hz] e frequenza centrale %.2f [Hz]', filter_type, B, f0);
title(titoloGrafico,'FontSize',16);
hold off;

% --- Calcolo della trasformata di Fourier dell'ingresso --- %
lunghezzaFft = 2^nextpow2(length(x));
X = fft(x,lunghezzaFft)*tempoC;
X = [X(lunghezzaFft/2+1:lunghezzaFft) X(1:lunghezzaFft/2+1)];
frequenza = freqC*linspace(-0.5,0.5,lunghezzaFft+1);

% --- Calcolo della trasformata di Fourier dell'uscita --- %
%lunghezzaFft=2^nextpow2(length(w));
W = fft(w,lunghezzaFft)*tempoC;
W = [W(lunghezzaFft/2+1:lunghezzaFft) W(1:lunghezzaFft/2+1)];
%frequenza=freqC*linspace(-0.5,0.5,lunghezzaFft+1);

nexttile; % aggiunge una "piastrella" alla figura
plot(frequenza/1e3,abs(X), 'Color', 'cyan', 'LineWidth', 1.5);
hold on;
plot(frequenza/1e3,abs(W), 'Color', 'black', 'LineWidth', 1.5);
grid on;
xlabel('Frequenza (kHz)', FontSize=14);
ylabel('Spettro di ampiezza', FontSize=14);
legend('|X(f)|', '|W(f)|', FontSize=12);
axis([0 8 0 1.2*max(abs(W))]);
title('Trasformate di Fourier dei segnali convoluti', 'FontSize',16)
hold off;

file = split(filename,".");
name = string(file(1));
out_filename = sprintf('./output/%s_output_%s_b%.0f_f%.0f.wav',name, filter_type, B, f0);
audiowrite(out_filename,[w.'*0.99/max(abs(w)),w.'*0.99/max(abs(w))],freqC);

end

---

Clipping

Lo script clip_audio prende tutto ciò che viene fatto in filtro_audio ed aggiunge un blocco non lineare a cui il segnale (relativo al file audio) viene dato in pasto a seguito del filtraggio. Questo blocco ha un andamento lineare finché |y|<yM. Superato questo limite, il valore di z satura a yM o –yM.

Lo scopo del blocco è quello di tagliare il segnale in ingresso, aggiungendovi dei disturbi udibili nel file di output generato al termine del programma. I disturbi sono inversamente proporzionali al valore di yM. Le trasformazioni non lineari aggiungono inoltre dei contributi in frequenza e contrastano quindi il lavoro dei filtri.

Per semplicità, clip_audio permette di applicare solamente un filtro passa-basso.

Dato il seguente blocco:

% --- Blocco non lineare --- %
switch block_number
    case 1 % --- blocco #1 --- %
        yM=0.12;      % ampiezza della saturazione
        z=y;          % copia non distorta del segnale di ingresso
        z(y>yM)=yM;   % clipping dei valori maggiori di yM
        z(y<-yM)=-yM; % clipping dei valori minori di -yM
    case 2
        yM=0.10;      % ampiezza della saturazione
        z=y;          % copia non distorta del segnale di ingresso
        z(y>0)=yM;    % clipping dei valori maggiori di 0
        z(y<0)=-yM;   % clipping dei valori minori di 0
    case 3
        yM=0.10;      % ampiezza della saturazione
        z=y;          % copia non distorta del segnale di ingresso
        z(y>0)=yM*(1-exp(-y(y>0)./yM)); % clipping dei valori maggiori di 0
        z(y<0)=-yM*(1-exp(y(y<0)./yM));% clipping dei valori minori di 0
    otherwise
        warning('Unexpected block number.')
        return;
end

Fornita in ingresso una sinusoide di ampiezza 0.2 e periodo uguale alla dimensione della banda, è interessante notare il confronto del segnale prima e dopo l’applicazione del blocco non lineare.

Applicando il 3° blocco distorsivo all’introduzione della canzone Scar Tissue dei Red Hot Chili Peppers, si ottengono i seguenti grafici.

Si può chiaramente apprezzare il disturbo dovuto al blocco non-lineare. I grafici mostrano in modo chiaro la distorsione, ma ascoltare la traccia audio generata dal comando audiowrite è sicuramente il modo migliore per comprendere l’entità della distorsione.

%% Utilizzo:
% clip_audio('Scar_Tissue.mp3', 3, 1.0e3)
% filename: nome del file da leggere (path, URL, ...)
% block_number: numero del blocco non-lineare da applicare
% B: dimensione di banda [Hz]
function clip_audio(filename, block_number, B)

% --- Lettura campioni --- %
freqC = 44.1e3;        % [Hz]: frequenza di campionamento = 44.1 kHz
tempoC = 1/freqC;      % [s] tempo di campionamento
durataTransitorio=0.1; % [s]
durata = 10.0;         % [s]
numeroCampioni = durata * freqC; % numero totale di campioni
inizioCampioni = durataTransitorio*freqC;
[xstereo,~]=audioread(filename,[inizioCampioni+1,inizioCampioni+1+numeroCampioni]);
x = (xstereo(:,1))'; % solo canale sinistro
tempo=0:tempoC:durata;
% x = 0.2 .* sin(B.*tempo);

% --- Filtro passa-basso --- %
T = 50/B; % durata risposta impulsiva [s]
tempoFiltro=0:tempoC:T;
tempoNorm = tempoFiltro-T/2;
h=2*B*sinc(2*B*tempoNorm).*rectpuls(tempoNorm/T);

% --- Filtraggio --- %
y=conv(x,h)*tempoC; % uscita filtrata
y=y(length(h):length(y));
tempoY=tempo(1:length(y))+T/2;
tempoY_plot=tempoY-durataTransitorio;

figure;
set(gcf,'defaultaxesfontname','Courier New')
plot(tempo-durataTransitorio, x, 'Color', 'cyan', 'LineWidth', 1.5);
hold on;
plot(tempoY_plot, y, 'Color', 'black', 'LineWidth', 1.5);
grid on;
xlabel('Tempo (s)','FontSize',12);
ylabel('Segnali temporali','FontSize',12)
legend('x(t)', 'y(t)');
axis([0 0.4 -0.8001 0.8001]);

% --- Blocco non lineare --- %
switch block_number
    case 1 % --- blocco #1 --- %
        yM=0.12;      % ampiezza della saturazione
        z=y;          % copia non distorta del segnale di ingresso
        z(y>yM)=yM;   % clipping dei valori maggiori di yM
        z(y<-yM)=-yM; % clipping dei valori minori di -yM
    case 2 % --- blocco #2 --- %
        yM=0.10;      % ampiezza della saturazione
        z=y;          % copia non distorta del segnale di ingresso
        z(y>0)=yM;    % clipping dei valori maggiori di 0
        z(y<0)=-yM;   % clipping dei valori minori di 0
    case 3 % --- blocco #3 --- %
        yM=0.10;      % ampiezza della saturazione
        z=y;          % copia non distorta del segnale di ingresso
        z(y>0)=yM*(1-exp(-y(y>0)./yM)); % clipping dei valori maggiori di 0
        z(y<0)=-yM*(1-exp(y(y<0)./yM));% clipping dei valori minori di 0
    otherwise
        warning('Unexpected block number.')
        return;
end

figure;
plot(tempoY_plot, y, 'Color', 'cyan', 'LineWidth', 1.5); hold on;
plot(tempoY_plot, z, 'Color', 'black', 'LineWidth', 1.5); grid on;
xlabel('Tempo (s)','FontSize',12);
ylabel('Segnali temporali','FontSize',12);
legend('y(t)', 'z(t)');
axis([0 0.4 -0.2001 0.2001]);

% --- Calcolo della trasformata di Fourier dell'ingresso --- %
lunghezzaFft=2^nextpow2(length(y));
Y=fft(y,lunghezzaFft)*tempoC;
Y=[Y(lunghezzaFft/2+1:lunghezzaFft) Y(1:lunghezzaFft/2+1)];
frequenza=freqC*linspace(-0.5,0.5,lunghezzaFft+1);

% --- Calcolo della trasformata di Fourier dell'uscita --- %
%lunghezzaFft=2^nextpow2(length(z));
Z=fft(z,lunghezzaFft)*tempoC;
Z=[Z(lunghezzaFft/2+1:lunghezzaFft) Z(1:lunghezzaFft/2+1)];
%frequenza=freqC*linspace(-0.5,0.5,lunghezzaFft+1);

figure;
plot(frequenza/1e3,20*log10(abs(Y)./max(abs(Y))), 'Color', 'cyan', 'LineWidth', 1.5);
hold on;
plot(frequenza/1e3,20*log10(abs(Z)./max(abs(Y))), 'Color', 'black', 'LineWidth', 1.5);
plot(frequenza/1e3,20*log10(abs(Y)./max(abs(Y))), 'Color', 'cyan', 'LineWidth', 0.5);
grid on;
xlabel('Frequenza (kHz)','FontSize',12);
ylabel('Spettro di ampiezza (dB)','FontSize',12);
legend('|Y(f)|', '|Z(f)|');
axis([0 8 -120 0]);

% scrittura del file audio a cui è stato applicato il filtro
% N.B. se il filename è un path relativo e inizia con uno o più punti,
% questo modo per estrapolare il nome senza estensione non funziona.
% Dunque 'filename' deve essere un path assoluto, a meno che il file non si
% trovi nella cartella nella quale è situato lo script
file = split(filename,".");
name = string(file(1));
out_filter= sprintf('./output/%s_output_filter_b%.0f.wav',name,B);
out_dist = sprintf('./output/%s_output_dist_b%.0f.wav',name,B);

audiowrite(out_filter,[y.',y.'],freqC);
audiowrite(out_dist,[z.',z.'],freqC);

Filtri applicati ad immagini

MATLAB è in grado di leggere il contenuto di un immagine mediante il comando imread, ad esempio:

x=imread('unife', '.tif');
% oppure, con singolo parametro:
x=imread('unife.tif');

Ad x viene assegnata la matrice bidimensionale contenente il segnale di intensità luminosa dell’immagine unife.tif presente nel path corrente. Le componenti a bassa frequenza descrivono le caratteristiche generali dell’immagini (ai quali l’occhio umano è più sensibile). Le componenti ad alta frequenza descrivono invece i dettagli dell’immagine.

Il comando imgaussfilt applica un filtraggio passa-basso di tipo Gaussiano. Si ha la possibilità di passare il parametro opzionale deviazione standard $\sigma$ (sigma):

sigma = 20; % default: 0.5
y1 = imgaussfilt(x, sigma);

Si può applicare un filtro passa-alto sottraendo al segnale originale il segnale ricavato dal filtraggio passa-basso.

z1=x-y1;

Si può ulteriormente filtrare l’immagine sottraendo il segnale filtrato ad un quadrato totalmente bianco.

quadratoBianco=255*ones(size(x),'uint8');
y2=quadratoBianco-z1;

Lo script che segue applica questi concetti per generare e mostrare le immagini filtrate, a partire da un’immagine di partenza fornita come parametro, è il seguente:

%% Utilizzo:
% filtro_lp_hp_immagine("Unife.tif")
function filtro_lp_hp_immagine(filename)

% estrapolazione di nome ed estensione dal parametro in ingresso
file = split(filename,".");
name = file(1);
ext = file(2);

x=imread(name, ext); % lettura immagine

figure;
imshow(x); % effettua il grafico dell'immagine
title('Immagine originale');

%% filtraggio passa-basso di tipo Gaussiano
sigma = 20;
y1 = imgaussfilt(x, sigma);

figure
imshow(y1);
title('Immagine elaborata con il filtro passa-basso');

% genera un'immagine a partire dal vettore bidimensionale in ingresso
% il nome di tale immagine è una concatenazione di stringhe
imwrite(y1,strcat(name,'_lp.', ext),ext);

%% Filtraggio passa-alto
z1=x-y1;

figure
imshow(z1);
title('Immagine elaborata con il filtro passa-alto');
imwrite(z1,strcat(name,'_hp.', ext),ext);

% si può filtrare l'immagine ottenuta sottraendo il segnale filtrato ad un quadrato totalmente bianco
quadratoBianco=255*ones(size(x),'uint8');
y2=quadratoBianco-z1;

figure
imshow(y2);
title('passa-alto (bianco)');
imwrite(y2,strcat(name,'_white_hp.', ext),ext);

A partire dall’immagine del Dipartimento di Ingegneria di Unife, sono illustrati i risultati dello script precedente:

Applicazioni della funzione fspecial

Tra gli innumerevoli modi per applicare filtri alle immagini, è bene parlare della funzione fspecial, la quale permette di applicare un filtro avente una determinata risposta impulsiva h passata come parametro. Per applicare un filtro con fspecial bisogna fare uso di imfilter, ad esempio:

# H = fspecial('nome-filtro', altri-parametri);

H = fspecial('motion',20,45);
B = imfilter(I,H);
imshow(B);

Estrazione delle componenti RGB

È possibile applicare altre tipologie di filtraggio, tra le quali l’estrazione delle componenti RGB (red, green, blue).

function rgb_extractor(filename)
% estrapolazione di nome ed estensione dal parametro in ingresso
% N.B. se il filename è un path relativo e inizia con uno o più punti,
% questo modo per estrapolare il nome senza estensione non funziona.
% Dunque 'filename' deve essere un path assoluto, a meno che il file non si
% trovi nella cartella nella quale è situato lo script
file = split(filename,".");
name = file(1);
ext = file(2);

my_image=imread(filename);

% componenti rossa, verde e blu in termini di intensità luminosa
red_brigh=my_image(:,:,1);
green_brigh=my_image(:,:,2);
blue_brigh=my_image(:,:,3);
black=zeros(size(my_image,1),size(my_image,2));
white=my_image(:,:,1:1:1);

% componenti rossa, verde e blu in termini di colore
red_component=cat(3,red_brigh,black,black);
green_component=cat(3,black,green_brigh,black);
blue_component=cat(3,black,black,blue_brigh);
just_black=cat(3,black,black,black);
just_white=cat(3,white,white,white);

figure
imshow(my_image);

figure
imshow(red_component);
imwrite(red_component,strcat('./output/', name,'_red.', ext),ext);

figure
imshow(green_component);
imwrite(green_component,strcat('./output/', name,'_green.', ext),ext);

figure
imshow(blue_component);
imwrite(blue_component,strcat('./output/', name,'_blue.', ext),ext);

figure
imshow(just_black);
figure
imshow(just_white);
end

La tabella che segue illustra i risultati ottenuti con lo script appena illustrato: