Is a high-passed signal the same as a signal minus a low-passed signal?

Question

My question is, if I want to high-pass a signal, is it the same as low-passing a signal and subtracting it from the signal? Is it theoretically the same? Is it practically the same?

I have searched (both on google and dsp.stackexchange) and I find conflicting answers. I have been playing with a signal and here are the results. I can't make much sense of it. Here is the signal with the sampling frequency once every four seconds. I designed a digital low-pass filter with the transition band from 0.8 mHz to 1 mHz and filtered the signal. I then also designed a high-pass filter with the same transition band and filtered the signal. Here are the results.

This first picture shows the original signal in black and the low-passed signal in blue. They are almost on top of each other but not quite. The red curve is signal minus the high-passed signal which is right on top of the signal.

This second image is just the first one zoomed in to show what's happening. Here we see that clearly the two are not the same. My question is why? Is it something about how I have implemented the two filters or is it something theoretical independent of my implementation? I don't know a lot about filter designing but I do know that it is notoriously counter-intuitive. Here is the full MATLAB code to reproduce all of this. I am using the filtfilt command to eliminate phase delays. But another thing to point out here is that the filters are not normalized. When I do sum(Hd.Numerator), I get 0.9930 for the low-pass and 0.007 for the high-pass. I don't see how to account for this. Should the output be scaled somehow because the coefficients don't add up to one? Could this scaling have something to do with this?

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')

Answer 1

In general you can't simply subtract a low-pass filtered version of a signal from the original one to obtain a high-pass filtered signal. The reason is as follows. What you're actually doing is implement a system with frequency response

$$H(\omega)=1-H_{LP}(\omega)\tag{1}$$

where $H_{LP}(\omega)$ is the frequency response of the low-pass filter. Note that $H_{LP}(\omega)$ is a complex function. What you probably want is

$$|H(\omega)|=\big|1-|H_{LP}(\omega)|\big|\tag{2}$$

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

Now write $H_{LP}(\omega)$ as

$$H_{LP}(\omega)=|H_{LP}(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase response of the low pass filter. If you defined a filter as

$$H_{HP}(\omega)=e^{j\phi(\omega)}-H_{LP}(\omega)=e^{j\phi(\omega)}\left(1-|H_{LP}(\omega)|\right)\tag{3}$$

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 $e^{j\phi(\omega)}$ 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

$$e^{j\phi(\omega)}=e^{-j\omega\tau}\tag{4}$$

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 $\tau$ in (4) an integer, which is achieved by choosing an even filter order $n$: $\tau=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 $x_{LP}[n]$
• delay the original signal by $\tau=n/2$ samples; let's call the delayed signal $x_d[n]$
• generate the high-pass filtered signal by subtraction: $x_{HP}[n]=x_d[n]-x_{LP}[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])

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])

Answer 2

Regarding scaling:

If you sum up the coefficients, you get the magnitude for DC. Hence, it makes sense that you get those numbers ($\approx 1$ for LP, $\approx 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.,

$|H_{LP}|+|H_{HP}| = 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.,

$|H_{LP}|^2+|H_{HP}|^2 = 1$