# Numerical error or mistake when filtering using Fourier?

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]


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]


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]


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]


• 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. – Chris Nadovich Jul 4 '17 at 3:36
• 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. – Chris Nadovich Jul 4 '17 at 3:38
• @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. – Hugh Jul 4 '17 at 8:46