7
$\begingroup$

First I define a function, just a sum of a few sin waves at different angular frequencies:

ubdat = 50;
ws = 10*{2, 5, 10, 20, 40}
fn = Table[Sum[Sin[w*x], {w, ws}], {x, 0, ubdat, .001}];
pts = Length@fn
ListPlot[fn, Joined -> True, PlotRange -> {{0, 1000}, All}]

{20, 50, 100, 200, 400}

enter image description here

If I take the Fourier transform and scale it correctly, you can see the correct peaks:

fnft = Abs@Fourier@fn;
fnftnormed = Table[{2*Pi*i/ubdat, fnft[[i]]}, {i, Length@fnft}];
ListPlot[fnftnormed, Joined -> True, PlotRange -> {{0, 500}, All}]

enter image description here

Now, I want to do a low pass filter on it, for, say, $\omega_c=140$. This should get rid of the peaks at 200 and 400, ideally. Doing it this way returns the same plots as above:

fnfilt = LowpassFilter[fn, 140];
ListPlot[fnfilt, Joined -> True, PlotRange -> {{0, 1000}, All}]
fnfiltft = Abs@Fourier@fnfilt;
fnfiltftnormed = 
  Table[{2*Pi*i/ubdat, fnfiltft[[i]]}, {i, Length@fnfiltft}];
ListPlot[fnfiltftnormed, Joined -> True, PlotRange -> {{0, 500}, All}]

I assume the problem is something to do with defining SampleRate, but the documentation explaining how it's defined or how to use it on the LowpassFilter page is very sparse:

By default, SampleRate->1 is assumed for images as well as data. For a sampled sound object of sample rate of r, SampleRate->r is used. With SampleRate->r, the cutoff frequency should be between 0 and $r*\pi$.

It appears to have a broken link at the bottom, so maybe that had something helpful. The page for SampleRate itself has even less info.

My naive attempt at choosing a sample rate would be dividing the number of samples by the total range, so in this case, Floor[pts/ubdat]=1000. Using this does affect the FT, but not a whole lot:

fnfilt = LowpassFilter[fn, 140, SampleRate -> 1000];
ListPlot[fnfilt, Joined -> True, PlotRange -> {{0, 1000}, All}]
fnfiltft = Abs@Fourier@fnfilt;
fnfiltftnormed = 
  Table[{2*Pi*i/ubdat, fnfiltft[[i]]}, {i, Length@fnfiltft}];
ListPlot[fnfiltftnormed, Joined -> True, PlotRange -> {{0, 500}, All}]

enter image description here

So what am I missing? I've tried googling for some sort of guide on using filters in Mathematica, but I can't find anything and it's very frustrating.

$\endgroup$
1
  • $\begingroup$ LowpassFilter is a black box, the transfer function is not known. The well-known and usual filters such as BiquadraticFilter, ButterworthFilter, ChebyshevFilter ... , the transfer function is known and you can build a filter based on the transfer function. $\endgroup$
    – user36273
    Jul 4, 2016 at 9:59

5 Answers 5

5
+50
$\begingroup$

I could not make LowpassFilter work for one-dimensional time histories either. I suggest you use a Butterworth filter which is more standard. Also I like to use random for testing because you can look at the transfer function. For the test below I add a few zeros at the end of the random to allow for filter ringing.

Here is your code modified for random.

ubdat = 50;
fn = Join[Table[RandomReal[{-1, 1}], {x, 0, ubdat - 10, .001}], 
   ConstantArray[0, Round[10/0.001]]];
pts = Length@fn
ListPlot[fn, Joined -> True, PlotRange -> {{0, 1000}, All}]

Mathematica graphics

The Fourier transform is less clear than with sine waves but stay with it

fnft = Abs@Fourier@fn;
fnftnormed = Table[{2*Pi*i/ubdat, fnft[[i]]}, {i, Length@fnft}];
ListPlot[fnftnormed, Joined -> True, PlotRange -> {{0, 500}, All}]

Mathematica graphics

Now I make a low pass filter, convert it to a time model and then use it on the data.

lowpass = ButterworthFilterModel[{3, 140}];
filt = ToDiscreteTimeModel[lowpass, 0.001];
fnfilt = RecurrenceFilter[filt, fn];
fnfiltft = Abs@Fourier@fnfilt;

fnfiltftnormed = 
  Table[{2*Pi*i/ubdat, fnfiltft[[i]]}, {i, Length@fnfiltft}];
ListPlot[fnfiltftnormed, Joined -> True, PlotRange -> {{0, 500}, All}]

Mathematica graphics

This is looking better with a clear reduction in level. Now we look at the relationship between the output and the input. I use your frequency axis.

ListLinePlot[
 Transpose[{fnfiltftnormed[[All, 1]], Abs[fnfiltft/fnft]}][[1 ;; 
    6000]], PlotRange -> All]

Mathematica graphics

Now there is a nice transfer function with a clear drop at 140 radians/sec.

If you repeat this with your sine waves you will get the result you want but you can't do the transfer function because the spectra are zero between the peaks.

Hope that helps.

$\endgroup$
6
  • $\begingroup$ Thank you Hugh, that's an interesting workaround. I'll try it if I can't get LPF[] to work, but I'd like to be able to use that. I imagine I must be missing something because filtering 1D data is pretty much the simplest application there is. My end goal is actually using BandstopFilter[] but I tried that as well and it didn't work. Thanks! $\endgroup$ Jun 28, 2016 at 15:00
  • $\begingroup$ I tried using random with Lowpassfilter and getting the transfer function but it did not work. This may be a bug but I suspect that this is intended for image work. There is a version of Butterworth that is a bandstop filter. I would suggest you use Butterworth. $\endgroup$
    – Hugh
    Jun 28, 2016 at 15:05
  • $\begingroup$ I guess I might have to. The reason I'm a little hesitant is that I'm unfamiliar with the workings of filters (beyond their basic concept), so I'll have to put in some time to figure out how it works. I notice that in your example, the rolloff isn't that sharp -- would it be a lot sharper if you just used a higher order than 3? $\endgroup$ Jun 28, 2016 at 15:09
  • $\begingroup$ Yes you can increase the order. Too high an order can become unstable. This is the way with filters the sharper the rolloff the worse the phase response and also the possibility of instability. $\endgroup$
    – Hugh
    Jun 28, 2016 at 15:12
  • $\begingroup$ @YungHummmma You may also be interested in a previous answer from Hugh himself on the topic of whether these are image processing or signal processing filters. (+1) $\endgroup$
    – MarcoB
    Jun 30, 2016 at 0:22
2
$\begingroup$

Actually, LowpassFilter is pretty easy to use, if you understand the parameters it needs. For instance, replacing the OPs line with

fnfilt = LowpassFilter[fn, 140, 200, SampleRate -> 1000];

gives exactly what you would expect.

ListPlot[fnfilt, Joined -> True, PlotRange -> {{0, 1000}, All}]
fnfiltft = Abs@Fourier@fnfilt;
fnfiltftnormed = 
  Table[{2*Pi*i/ubdat, fnfiltft[[i]]}, {i, Length@fnfiltft}];
ListPlot[fnfiltftnormed, Joined -> True, PlotRange -> {{0, 500}, All}]

enter image description here

LowpassFilter works fine in 1D as well as 2D. The filter needs to know the cutoff frequency (140 in this case), the sampling rate (1000 in the OPs example), and it needs help choosing the length of the filter. Longer filter lengths allow sharper cutoffs but are computationally more expensive.

$\endgroup$
2
  • $\begingroup$ Thanks, I like this simple solution. But I still don't really know what SampleRate does here, or how to choose it, which makes me nervous. For example, using your solution, if I choose SampleRate->2000 instead, it starts cutting into frequencies below the cutoff frequency. Why would increasing the sample rate do this? $\endgroup$ Jul 7, 2016 at 15:26
  • $\begingroup$ The sample rate is the number of samples per second. You generated your signal using fn = Table[Sum[Sin[w*x], {w, ws}], {x, 0, ubdat, .001}];. Thus you have chosen to take .001 samples each unit of x, hence you have already defined the sample rate, you just need to let Lowpass know. $\endgroup$
    – bill s
    Jul 7, 2016 at 15:29
1
$\begingroup$

As MarcoB remembered I have looked at this before and concluded that LowpassFilter is for image processing not for signal processing. Lets go through your example with random noise and see what we get. I start by making some random noise with zeros before and after.

sr = Round[1/0.001];
fn = Join[ConstantArray[0, 500], RandomReal[{-1, 1}, 49000], 
   ConstantArray[0, 500]];
pts = Length@fn
ListPlot[fn, Joined -> True, PlotRange -> {{0, 1000}, All}]

Mathematica graphics

Now do the filtering. Plot the start of the original and the filtered time history against point number.

fnfilt = LowpassFilter[fn, 140, SampleRate -> sr];
ListLinePlot[{fn, fnfilt}, PlotRange -> {{475, 525}, All}]

Mathematica graphics

The filtered signal begins before the original signal. This is a non-causal response which is usually avoided in signal processing. At the end there is also a signal after the original signal has stopped (not plotted). This is acceptable in signal processing and is due to ringing of the filter.

Take the Fourier transforms and then plot the filtered spectrum and the ratio of the filtered spectrum to the original spectrum. This is the filter transfer function. Also make the frequency axis. Note that the first point is zero frequency.

fnft = Abs@Fourier@fn;
fnfiltft = Abs@Fourier@fnfilt;
freqs = Table[2.*Pi*(i - 1) sr/(pts - 1), {i, pts}];
ListLinePlot[Transpose[{freqs, fnfiltft}], 
 PlotRange -> {{0, 500}, All}]
ListLinePlot[Transpose[{freqs, fnfiltft/fnft}], 
 PlotRange -> {{0, 500}, All}, 
 Epilog -> {Pink, Line[{{0, 1/Sqrt[2]}, {140, 1/Sqrt[2]}, {140, 0}}]}]

Mathematica graphics

Mathematica graphics

I have added pink lines to show the standard cut off-frequency. As a filter transfer function goes this is not very good. It is not flat in the pass band and does not drop to the usual 1/Sqrt[2] at the cut-off frequency.

I really recommend that you use the signal processing software I described here

$\endgroup$
1
$\begingroup$

I've found that specifying the kernel filter length in LowpassFilter will give the better performance than leaving it unspecified. I'll show what I came up with. Note that I've modified your problem somewhat. I use Fourier with the FourierParameters set for "signal processing". I renamed some of your variables for clarity. I also shortened the length of the sampled signal to 4096 samples. This shouldn't affect anything but the resolution of the frequency bins that result from calling Fourier. Then, I moved the frequency of one of the Sin functions from 100 rad/sec to 140 rad/sec. This makes it easier to see the effect of the filter at the cutoff frequency. All of the DFT frequency response plots use dB with the normalize frequency data. This makes it easier to see the relative drop off of the filter.

Here's my redefinition for Fourier:

myFourier[v_ /; VectorQ[v] || MatrixQ[v]] := Fourier[v, FourierParameters -> {1, -1}]

This generates the unfiltered signal and sets up the others variables:

signalDataLength = 4096;
lpfCutoffAngular = 140;
wSignal = 10*{2, 5, 14, 20, 40};
sampleRate = 1000.0;
sampleRateAngular = 2 Pi sampleRate;
unfilteredSignal = Table[Sum[Sin[w*(x/sampleRate)], {w, wSignal}], {x, 0,signalDataLength - 1}];
ListLinePlot[unfilteredSignal, PlotRange -> {{0, 512}, All},AxesLabel -> "Unfiltered Signal"]

Unfiltered Signal

Here's the normalized DFT of the unfiltered signal:

unfilteredDFT = Abs[myFourier[unfilteredSignal]];
angularFreqAxis = sampleRateAngular Range[0, signalDataLength/2 - 1]/signalDataLength;
unfilteredData = Transpose@{angularFreqAxis, 20 Log10[unfilteredDFT[[;; signalDataLength/2]]/Max[unfilteredDFT]]};
ListLinePlot[unfilteredData, PlotRange -> {{0, sampleRateAngular/10}, {-50, 0}}, AxesLabel -> {"rad/sec", "dB"}]

Unfiltered DFT

Now, select the filter kernel length. I ended up using 256 since it gave a the filter a good match to your requirements of a cutoff of 140 rad/sec:

filterKernelLength = 256;
lpfSignal = LowpassFilter[unfilteredSignal, lpfCutoffAngular, filterKernelLength, SampleRate -> Round[sampleRate]];
lpfSignalDFT = Abs[myFourier[lpfSignal]];
lpfSignalData = Transpose@{angularFreqAxis, 20 Log10[lpfSignalDFT[[;; signalDataLength/2]]/Max[lpfSignalDFT]]};

ListLinePlot[{unfilteredData, lpfSignalData}, PlotRange -> {{0, sampleRateAngular/10}, {-50, 0}}, Joined -> True, Epilog -> {Red, Line[{{0, -3}, {lpfCutoffAngular, -3}, {lpfCutoffAngular, -100}}]}, PlotLegends -> {"Unfiltered", "Low Pass Filtered"}, AxesLabel -> {"rad/sec", "dB"}]

Unfiltered and Filtered DFT's

You should be able to see that the filtered response at 200 and 400 rad/sec is significantly attenuated.

Here's a plot of a portion of the unfiltered and filtered signals:

ListLinePlot[{unfilteredSignal, lpfSignal}, PlotRange -> {{0, 512}, ll}, PlotLegends -> {"Unfiltered", "Low Pass Filtered"}]

enter image description here

To look at the filter's impulse response, I did the following:

impulseData = SparseArray[{{signalDataLength/2} -> 1}, {signalDataLength}];
impResTS = TimeSeries[ impulseData, { Range[0, signalDataLength - 1]/sampleRate}];
impulseResponse = LowpassFilter[impResTS, Quantity[140, "Radians")/("Seconds")], filterKernelLength];

ListLinePlot[impulseResponse, PlotRange -> All, PlotLabel -> "Low Pass Filter Impulse Response"]

Low Pass Filter Impulse Response

impulseResponseDFT = Abs[myFourier[impulseResponse["Values"]]];
impulseResponseDFTData = Transpose@{angularFreqAxis, 20 Log10[impulseResponseDFT[[;; signalDataLength/2]]/Max[impulseResponseDFT]]};

ListLinePlot[tData, Joined -> True, PlotRange -> {{0, sampleRateAngular/10}, {-80, 0}}, Epilog -> {Red, Line[{{0, -3}, {lpfCutoffAngular, -3}, {lpfCutoffAngular, -100}}]}, AxesLabel -> {"rad/sec", "dB"}, PlotLabel -> "Low Pass Filter Frequency Response"]

The filter cutoff is slightly more the 3 dB down at 140 rad/sec. That should meet your requirements.

Low Pass Filter Impulse Frequency Response

$\endgroup$
1
$\begingroup$

I choose a simple Butterworth filter and show, based on the spectra, how it works. Cutoff frequency and filter order can be set.

ubdat = 50;
ws = 10*{2, 5, 10, 20, 40};
fn = Table[Sum[Sin[w*x], {w, ws}], {x, 0, ubdat, 0.001}];

Manipulate[
 filter = ButterworthFilterModel[{"Lowpass", filterOrder, cutoffFreqency}];
 discreteFilter = ToDiscreteTimeModel[filter, 0.001];

(* for 1-D filter one can set `OutputResponse` instead of `RecurrenceFilter` *)
(* recurrence = RecurrenceFilter[discreteFilter, fn]; *)
recurrence =OutputResponse[discreteFilter, fn][[1]]; 

 fourier = Abs@Fourier@recurrence; 
 fourierNormed = Table[{2*Pi*i/ubdat, fourier[[i]]}, {i, Length@fourier}];

 ListLinePlot[fourierNormed, PlotRange -> {{0, 500}, All}],
 {filterOrder, {1, 2, 3, 4, 5}},
 {cutoffFreqency, {400, 300, 200, 140, 100}}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.