Czy sygnał górnoprzepustowy jest taki sam jak sygnał minus sygnał dolnoprzepustowy?

14

Moje pytanie brzmi: jeśli chcę przekazać sygnał górnoprzepustowo, czy to to samo, co przekazanie sygnału dolnoprzepustowego i odjęcie go od sygnału? Czy to teoretycznie to samo? Czy to praktycznie to samo?

Szukałem (zarówno w Google, jak i dsp.stackexchange) i znajduję sprzeczne odpowiedzi. Bawiłem się sygnałem i oto wyniki. Nie mam większego sensu. Oto sygnał z częstotliwością próbkowania raz na cztery sekundy. Zaprojektowałem cyfrowy filtr dolnoprzepustowy z pasmem przejściowym od 0,8 MHz do 1 MHz i przefiltrowałem sygnał. Następnie zaprojektowałem również filtr górnoprzepustowy z tym samym pasmem przejściowym i przefiltrowałem sygnał. Oto wyniki.

enter image description here

To pierwsze zdjęcie pokazuje oryginalny sygnał w kolorze czarnym, a dolnoprzepustowy sygnał w kolorze niebieskim. Są prawie na sobie, ale nie całkiem. Czerwona krzywa to sygnał minus górnoprzepustowy sygnał, który znajduje się bezpośrednio nad sygnałem.

enter image description here

To drugie zdjęcie jest tylko pierwszym powiększonym obrazem pokazującym, co się dzieje. Tutaj widzimy, że wyraźnie te dwa nie są takie same. Moje pytanie brzmi: dlaczego? Czy jest to związane z tym, jak zaimplementowałem dwa filtry, czy jest to coś teoretycznego niezależnego od mojej implementacji? Nie wiem dużo o projektowaniu filtrów, ale wiem, że jest to sprzeczne z intuicją. Oto pełny kod MATLAB do odtworzenia tego wszystkiego. Używam polecenia filtfilt, aby wyeliminować opóźnienia fazowe. Ale inną rzeczą, na którą należy tutaj zwrócić uwagę, jest to, że filtry nie są znormalizowane. Kiedy robię sumę (Hd.Numerator), dostaję 0,9930 dla dolnego przejścia i 0,007 dla górnego przejścia. Nie wiem, jak to wytłumaczyć. Czy wynik powinien być jakoś skalowany, ponieważ współczynniki się nie sumują? Czy to skalowanie może mieć z tym coś wspólnego?

close all
clear all
clc

data = dlmread('data.txt');

Fs    = 0.25;    % Sampling Frequency
N     = 2674;    % Order
Fpass = 0.8/1000;  % Passband Frequency
Fstop = 1/1000;   % Stopband Frequency
Wpass = 1;       % Passband Weight
Wstop = 1;       % Stopband Weight
dens  = 20;      % Density Factor

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop], {dens});
Hd = dsp.FIRFilter('Numerator', b);
sum(Hd.Numerator)
datalowpassed = filtfilt(Hd.Numerator,1,data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs    = 0.25;    % Sampling Frequency
N     = 2674;    % Order
Fstop = 0.8/1000;  % Stopband Frequency
Fpass = 1/1000;   % Passband Frequency
Wstop = 1;       % Stopband Weight
Wpass = 1;       % Passband Weight
dens  = 20;      % Density Factor

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, [0 Fstop Fpass Fs/2]/(Fs/2), [0 0 1 1], [Wstop Wpass], {dens});
Hd = dsp.FIRFilter('Numerator', b);
sum(Hd.Numerator)
datahighpassed = filtfilt(Hd.Numerator,1,data);

figure
subplot(2,1,1)
plot(data,'-ko')
hold on
plot(datalowpassed,'-bo')
plot(data-datahighpassed,'-ro')
legend('Original Signal','Low-Passed','Signal - High-Passed')
subplot(2,1,2)
plot(data-datalowpassed,'-bo')
hold on
plot(datahighpassed,'-ro')
legend('Signal - Low-Passed','High-Passed')
Fixed Point
źródło
1
Twoja kolejność filtrów jest bardzo wysoka. Możliwe, że napotkasz problemy numeryczne podczas procesu projektowania. Sprawdź zaprojektowane filtry, obliczając i wykreślając wielkość FFT. Ponadto zobacz moją odpowiedź poniżej, aby dowiedzieć się, jak używać odejmowania do generowania filtra górnoprzepustowego z filtra dolnoprzepustowego.
Matt L.
W ten sposób można wygenerować sygnał górnoprzepustowy, ale kolejność wycofywania jest zawsze pierwsza: sound.westhost.com/articles/derived-xovers.htm Jeśli najpierw odłożysz sygnał o opóźnienie grupowe LPF przed odjęciem , możesz dostać 3. zamówienie
endolith

Odpowiedzi:

8

Zasadniczo nie można po prostu odjąć filtrowanej dolnoprzepustowej wersji sygnału od oryginalnej, aby uzyskać filtrowany dolnoprzepustowy sygnał. Powód tego jest następujący. To, co faktycznie robisz, to wdrożenie systemu z pasmem przenoszenia

(1)H(ω)=1HLP(ω)

where HLP(ω) is the frequency response of the low-pass filter. Note that HLP(ω) is a complex function. What you probably want is

(2)|H(ω)|=|1|HLP(ω)||

but this is generally not the case when (1) is satisfied.

Now write HLP(ω) as

HLP(ω)=|HLP(ω)|ejϕ(ω)

where ϕ(ω) is the phase response of the low pass filter. If you defined a filter as

(3)HHP(ω)=ejϕ(ω)HLP(ω)=ejϕ(ω)(1|HLP(ω)|)

then this new filter would indeed satisfy the magnitude relation (2). So if you do the subtraction as in (3), you can realize a high-pass filter with a magnitude response that is complementary to the response of the low pass filter. This means that you need to filter the signal by an all-pass filter with frequency response ejϕ(ω) before subtracting from it the low-pass filtered version of the signal.

In practice this is very simple if the low-pass filter has a linear phase response, because then the phase term is given by

(4)ejϕ(ω)=ejωτ

and corresponds to a simple delay. FIR filters can have an exactly linear phase response, and the ones designed by the Parks-McClellan algorithm do have a linear phase. It is advantageous to make the delay τ in (4) an integer, which is achieved by choosing an even filter order n: τ=n/2

So what you have to do is the following:

  • design a linear phase FIR low-pass filter with an even order
  • filter the signal using filter(); let's call the result xLP[n]
  • delay the original signal by τ=n/2 samples; let's call the delayed signal xd[n]
  • generate the high-pass filtered signal by subtraction: xHP[n]=xd[n]xLP[n]

Here is a very simple illustration in Matlab/Octave

h_lp = fir1(100,.3);                        % low-pass design
h_hp = [zeros(50,1);1;zeros(50,1)] - h_lp;  % high-pass design by subtraction
[H_lp,w] = freqz(h_lp,1,1024);
[H_hp,w] = freqz(h_hp,1,1024);
plot(w/2/pi,20*log10(abs(H_lp)),w/2/pi,20*log10(abs(H_hp)))
grid, axis([0,.5,-100,5])

enter image description here

EDIT:

Since you used the filtfilt command, the phase is artificially eliminated, and Equations (1) and (2) above become equivalent, because all frequency responses are in fact squared magnitudes of the designed responses. So, apart from small differences due to the filter design process, numerical errors, and small differences caused by the filtfilt function (automatically chosen initial conditions minimizing start and end transients), the result of subtracting filtered data from the original data should closely resemble direct filtering with a complementary filter. Since this is not the case in your example, I suspect that the filter design routine gives you trouble due the extremely high filter order. If you do the same with simpler filters (I chose n=100), you get what you would expect. In the figure below you see a section of the data in blue, the output of the low pass filter in green, and the result of the subtraction of the output of the high pass filter from the original data in red. The green and red curves are virtually identical.

x = load('data.txt');           % data to be filtered
h_lp = fir1(100,.3);            % LP impulse response
h_hp = fir1(100,.3,'high');     % HP impulse response
y = filtfilt(h_lp,1,x);         % apply low pass filter
yh = filtfilt(h_hp,1,x);        % apply high pass filter
yd = x - yh;                    % low pass by difference with high pass filter
n = 1:length(x);
plot(n,x,n,y,'g.',n,yd,'r')
axis([3500,4000,140,150])

enter image description here

Matt L.
źródło
If you are trying to design a high pass filter this way, you have to be careful with the specifications of the low-pass filter. The stop band attenuation in the low-pass filter is usually high enough to achieve small pass band ripple in the high pass, but the pass band ripple in the LP filter often does not achieve sufficient stop band attenuation in the HP filter.
David
Thanks for the detailed response. It cleared up a few things.
Fixed Point
3

Regarding scaling:

If you sum up the coefficients, you get the magnitude for DC. Hence, it makes sense that you get those numbers (1 for LP, 0 for HP).

In addition to Matt L.'s excellent answer one can just point out that what he is using is referred to as magnitude complementary filters, which is the common case for linear-phase FIR filters, i.e.,

|HLP|+|HHP|=1

When creating filters from two parallel allpass sections and adding/subtracting the outputs, the lowpass/highpass filters will instead be power complementary, i.e.,

|HLP|2+|HHP|2=1

Oscar
źródło