0
$\begingroup$

I am comparing the filtering method using a recurrence filter with filtering using Fourier. I think both methods should give the same answer but they differ by a small amount that could be numerical noise but seems to be too large. Am I making a mistake or could the noise be this large? Here is a simple example where I filter a sine wave. I add plenty of zeros at the end so that any ringing in the filter is properly accounted for.

sr = 5000.; (* sample rate *)
t0 = -0.2 + 1/sr;   (* time before start, ordinates = 0 *)
t1 = 0.3; (* end of signal *)
t2 = 4.;  (* time after signal, ordinates = 0 *)
f0 = 123.4; (* frequency of sine wave *)
f1 = 110; (* start filter frequency *)
f2 = 140;  (* end filter frequency *)
a = Join[Transpose[{Range[t0, 0 - 1/sr, 1/sr], 
    ConstantArray[0, Round[-t0 sr]]}], 
  Table[{t, Sin[2 π f0 t]}, {t, 0, t1, 1/sr}], 
  Transpose[{Range[t1 + 1/sr, t2, 1/sr], 
    ConstantArray[0, Round[(t2 - t1) sr]]}]];
ListLinePlot[a, PlotRange -> All, Frame -> True, AspectRatio -> 1/4]

Mathematica graphics

Now I make the filter and then filter the input. This all looks good.

filt = ToDiscreteTimeModel[
   ButterworthFilterModel[{"Bandpass", 4, 2 \[Pi] {f1, f2}}], 1/sr];
b = RecurrenceFilter[filt, a[[All, 2]]];
ListLinePlot[b, PlotRange -> All, Frame -> True, AspectRatio -> 1/4]

Mathematica graphics

Now we do the Fourier method. I start by making the frequency response function from the filter. Then I multiply this by the Fourier transform of the input and then finally I take the inverse Fourier transform of the product. The result is similar but the error is large. The result also has a small imaginary part which should be zero.

ClearAll[f];
n = Length[a];
{{h}} = filt[E^(I 2 π f/sr)];
frf = Table[h, {f, 0, sr - sr/n, sr/n}];
x = Fourier[a[[All, 2]], FourierParameters -> {-1, -1}];
y = InverseFourier[frf x, FourierParameters -> {-1, -1}];
ListLinePlot[Re@y, PlotRange -> All, Frame -> True, AspectRatio -> 1/4]
ListLinePlot[Re@(b - y), PlotRange -> All, Frame -> True, 
 AspectRatio -> 1/4]
ListLinePlot[Im@y, PlotRange -> All, Frame -> True, 
 AspectRatio -> 1/4]

Mathematica graphics Mathematica graphics Mathematica graphics

I would not expect numerical errors of the order of 10^-4 I would expect them to be the order of (10^-10). I suspect the Fourier method rather than the filter method because of the appearance of an imaginary part. Is there a mistake or have numerical errors crept in?

Edit

I have just discovered that I can get the errors down by one more order of magnitude if I use cos + I sin rather than an exponential function. This continues to point towards a numerical problem. How can I improve the numerics further?

ClearAll[f];
n = Length[a];
{{h}} = filt[Cos[2 π f/sr] + I Sin[2 π f/sr]];
frf = Table[h, {f, 0, sr - sr/n, sr/n}];
x = Fourier[a[[All, 2]], FourierParameters -> {-1, -1}];
y = InverseFourier[frf x, FourierParameters -> {-1, -1}];
ListLinePlot[Re@y, PlotRange -> All, Frame -> True, AspectRatio -> 1/4]
ListLinePlot[Re@(b - y), PlotRange -> All, Frame -> True, 
 AspectRatio -> 1/4]
ListLinePlot[Im@y, PlotRange -> All, Frame -> True, 
 AspectRatio -> 1/4]

Mathematica graphics Mathematica graphics Mathematica graphics

$\endgroup$
  • $\begingroup$ If your frequency domain functions have exactly even real parts and exactly odd imaginary parts, their transform should be purely real. You say there's a small imaginary part. How small? Are they really exactly even and odd. You should check. If you're off by one it'll mess everything up. $\endgroup$ – Chris Nadovich Jul 4 '17 at 3:36
  • $\begingroup$ Also, recursive filters like Butterworth can be noisy, and can limit cycle (output with no input). Have you thought about this aspect? Recursive filters are rarely used for this very reason. Most people use FIR. $\endgroup$ – Chris Nadovich Jul 4 '17 at 3:38
  • $\begingroup$ @ChrisNadovich Thanks for these comments I just checked x which seems fine but the frf may be off. I will report back. This could be a good hint for the problem. $\endgroup$ – Hugh Jul 4 '17 at 8:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.