A Simple Lowpass Filter Behaving so Strangely

Question

I don't know much about filter design. After reading a few articles on wikipedia I am playing with filter design and I see something weird and would like some help in understanding what exactly is going on here. So I have a signal measured uniformly in time (every 30 seconds) with 1006 measurements in total and I want to estimate its power spectrum density eventually. However I want to look only at certain high frequencies so I decided to run it through a low pass filter to get the power in the zero channel and the neighboring small frequencies. Then I will subtract the filtered signal and look at the residual for power at the higher frequencies.

First thing, I just wrote my own naive ideal low pass filter where I take the FFT, zero out all frequencies above 0.8 mHz and then take the inverse FFT and this is the picture I get.

This picture makes perfect sense to me because the filtered signal just oscillates around the original signal and those are the low frequencies I see.

Next I used the GUI based filter builder included in MATLAB's DSP System Toolbox to design my own filter. So I made a low pass, finite impulse response, single-rate filter using the equiripple method and the minimum order possible with the passed frequency amplitude being 1dB and the stopped frequency amplitude attenuated 60dB. Since the range must start with a strictly positive number I used the allowed range of 10e-9 Hz to 0.0008 Hz. The filter is implemented as a direct-form FIR. Here is the frequency response.

My question is when I run the same data through this designed filter, the filtered signal looks so strange.

Can someone please enlighten me as to what is going on here? Since I am using MATLAB's own implementation I don't think there is anything wrong with the math or the computations so it must be my own understanding. Is this expected? Is it the fact that there is nothing wrong anywhere and the picture is simply what it is? Is there a different filter I should be looking at instead or something? I was expecting this pic to look more or less like the first one where some low frequencies with the background power are just overlayed.

I played with high pass filters and even band-pass filters over the specific frequency range I want to look at but they all behave strangely. Is it something to do with the signal itself? The background power is huge so is it leaking all over the spectrum and messing everything up or something? Thanks!

---

Edit 1

Thanks everyone for your comments. My original signal is that blue U shaped signal. I can post the measurements if you guys want here as soon as I figure out how to make a collapsable section. The original problem is to estimate the power spectral density between 0.8mHz and 8mHz. My first instinct was to just have a low pass filter and then subtract it. And then use the multitaper method on the residual to estimate the PSD. Subtraction came to mind because of the detrending. Since one usually detrends the data, I decided to detrend using the filtered signal.

First I wrote my own naive digital filter, filtering in the frequency space. Here's the MATLAB code.

Fs = 1/30;         % Sampling Frequency - per second
L = length(data);  % Length of the signal
fftb = fft(data);
f = Fs/2*linspace(0,1,L/2+1);
f(end+1:2*end-1)=fliplr(f(2:end));
fftb(f(:)>0.0008)=0;
naivefiltered = real(ifft(fftb));

Just take the FFT, zero the frequencies I don't want, then take the inverse FFT which sometimes is complex for some reason so I take the real part and I get the first picture up above.

Then I used MATLAB's filter designing tool to design this filter and here is the code for that.

Fpass = 3.3135e-09;  % Passband Frequency
Fstop = 0.0008;      % Stopband Frequency
Apass = 1;           % Passband Ripple (dB)
Astop = 60;          % Stopband Attenuation (dB)
Fs    = 1/30;        % Sampling Frequency
h = fdesign.lowpass('fp,fst,ap,ast', Fpass, Fstop, Apass, Astop, Fs);
Hd = design(h, 'equiripple','MinOrder', 'any','StopbandShape', 'flat');

Then for the third picture, I was using

smbtotal = filter(Hd,data);

but after endolith's comment (I didn't even know about filtfilt) I read about filtfilt and see that "filter()" introduces a huge time lag so when I use

smbtotal = filtfilt(Hd.Numerator,1,data);

the picture does look much better. Here it is.

So it looks like this time/phase shift was the problem. So now I have three questions.

1.For the type of data you see here (U shaped ranging over several orders of magnitudes) if I want to estimate the power at a specific band (0.8mHZ to 8mHz), which technique is better? Should I subtract a low-pass filtered signal? Should I high-pass filter the signal? Or should I band-pass filter it and then estimate the PSD?

2.Designing filters is non-trivial. But why is applying filters non-trivial? Why does using "filter()" introduce this time/phase shift?

3.Does it matter if I convolve in physical space or multiply in frequency domain? Which one is better? It looks to me like if I multiply in frequency domain (like I am doing with my naive filter) there is no issue with time/phase shifting. But if I convolve in physical space (like using "filter()" in MATLAB) then you have a problem with time/phase shifting? So is it always better to just multiply in frequency space?

Thanks everyone for your time. Appreciate it!

Answer 1

Real filters for real-time applications approximate the ideal filter by truncating and windowing the infinite impulse response to make a finite impulse response; applying that filter requires delaying the signal for a moderate period of time, allowing the computation to "see" a little bit into the future. This delay is manifested as phase shift.

You can calculate the phase shift manually(i.e by hand by following these steps).

1. Compute the Fourier transform of the filter.
2. Calculate the angle(or theta) of the Fourier transform. This means F(w) = |F(w)|(theta)
3. The theta is the phase shift of the signal.

You can then compute the time delay by multiplying the phase shift angle with the sampling(?) frequency.

If you want to work without time delays(or phase shifts) you should look at zero-phase filters which can be formed via bi-directional filtering.

When you're multiplying in frequency domain, you're ignoring certain aspects of the time-domain signal and therefore not obtaining the true representation of your filter.

For example: In the frequency domain I can easily define F(w) such that it is 0 for 0:Fcutoff and 1 for Fcutoff:infinity.

Sample algo for non-ideal High-pass filter in frequency domain would be :

 // Return RC high-pass filter output samples, given input samples,
 // time interval dt, and time constant RC
 function highpass(real[0..n] x, real dt, real RC)
   var real[0..n] y
   var real α := RC / (RC + dt)
   y[0] := x[0]
   for i from 1 to n
     y[i] := α * y[i-1] + α * (x[i] - x[i-1])
   return y

Therefore avoid designing the high-pass filters in the Fourier domain alone. Here alpha signifies the delay which you would have to artificially introduce.

This would be an ideal high-pass filter with magnitude 1, which is un-realizable in real life because of the complexities that sampling brings into the picture.

Answer 2

You mention that you want to interrogate the frequency response of your signal, specifically, that you want to know the levels of your signals high frequency components. (From: "However I want to look only at certain high frequencies").

In this case, if you just compute the DFT of your signal using MATLAB's built in function, this will tell you the frequency response of your signal. For example, you can do:

freqz(signal, 1, 2^nextpow2(length(signal)), your_sampling_frequency_here);

This should solve your problem as I understand it. From there, you can simply study the high-frequency components of your signal.

As a side note, you cannot, generally, filter a signal, then subtract the filtered version and look at the residual. This is because the very act of filtering your signal will introduce time/phase delays, which means that a subtraction between your filtered signal and original signal will be distorted accordingly. There are some techniques you could use to take those delays into account, but as a general case, be wary of doing such things.