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!
filtfilt()
? $\endgroup$