Ideal/Brick wall Low/High-pass filters: Fourier approach not working

Question

I am trying to take a signal (as a list) and separate it into the low frequency and high frequency parts around some split frequency f.

Given that what little signal processing I ever studied is now decades in the past, I'm probably doing the wrong thing anyway but I am puzzled by the results and maybe someone can explain them and how I should do it.

I thought the easy way would be: Fourier my signal, clear all the components above (below) f and InverseFourier back...

The test signal is a 1s chord of C (CEG in equal parts C @261Hz) for f=300 the low part seems to give the C, but the high part is not just E & G as the following graphs show.

The code (including debug for completeness) is as follows

oneWave[f_, sampleRate_] := 
 Table[Cos[i 2 Pi], {i, 0, (f - 1)/f, f/sampleRate}]
nWaves[f_, sampleRate_, n_] := 
 Flatten[Table[oneWave[f, sampleRate], {i, 1, n}]]
ceg4 = ListPlay[a = nWaves[261.63`, 44100, 262]; 
  1/3 (a + Take[nWaves[329.63`, 44100, 330], Length[a]] + 
     Take[nWaves[392.`, 44100, 392], Length[a]]), SampleRate -> 44100]

requencySplitDebug[signalList_List, splitF_?NumberQ, myDebugParams_List]:= 
    Module[
        (* "Some common choices for {a,b} are {0,1} (default), {-1,1} (data analysis), {1,-1} (signal processing). " *)
        {signalLen, lfPart, hfPart, fourierData, myDebugGraphics, myDebugData, lfIF, hfIF, fourierA=0, fourierB=1,rPad, lPad},
        signalLen = Length[signalList];
        fourierData = Fourier[signalList, FourierParameters->{fourierA, fourierB}];
        rPad = PadRight[Take[fourierData,{1,splitF}],signalLen];
        lPad = PadLeft[Take[fourierData,{splitF + 1,signalLen}],signalLen];
        lfIF = InverseFourier[rPad,FourierParameters->{fourierA, fourierB}];
        hfIF = InverseFourier[lPad,FourierParameters->{fourierA, fourierB}];
        lfPart = Chop[Re[lfIF]];
        hfPart = Chop[Re[hfIF]];
        If[Length[myDebugParams]==3 && myDebugParams[[1]]==True && myDebugParams[[3]] > myDebugParams[[2]],
            Print["Debugging On"];
            myDebugGraphics = 
                GraphicsGrid[{
                    {ListLinePlot[Take[signalList,{myDebugParams[[2]],myDebugParams[[3]]}], PlotRange->All,PlotLabel->"Signal Input Detail"],
                    ListLinePlot[Take[Abs[fourierData],{myDebugParams[[2]],myDebugParams[[3]]}], PlotRange->All, PlotLabel->"Signal Input Spectrum"]}, 
                    {ListLinePlot[Take[lfPart,{myDebugParams[[2]],myDebugParams[[3]]}], PlotRange->All,PlotLabel->"Lower Part Signal Detail"],
                    ListLinePlot[Take[Abs[Fourier[lfPart,FourierParameters->{fourierA, fourierB}]],{myDebugParams[[2]],myDebugParams[[3]]}], PlotRange->All, PlotLabel->"Lower Part Spectrum"]},
                    {ListLinePlot[Take[hfPart,{myDebugParams[[2]],myDebugParams[[3]]}], PlotRange->All,PlotLabel->"Upper Part Signal Detail"],
                    ListLinePlot[Take[Abs[Fourier[hfPart,FourierParameters->{fourierA, fourierB}]],{myDebugParams[[2]],myDebugParams[[3]]}], PlotRange->All, PlotLabel->"Upper Part Spectrum"]}
                }];
            myDebugData = {fourierData, lfIF, hfIF};
        ];

        If[Length[myDebugData]==0,
            Return[{lfPart,hfPart}],
            Return[{lfPart,hfPart,myDebugGraphics, myDebugData}]
        ];
    ];

asplit = frequencySplitDebug[AudioData[ceg4], 300, {True, 1, 600}];
Show[asplit[[3]]]

I really do want to get as close as possible to an ideal filter with execution time on a modest PC (i3) of ~1s for a list of length ~1000 because my real signal is data not audio and both parts need to be as clean as possible for the next processing steps.

I also tried a 7th order Buttwerworth filter on ceg4 (40k points) with high attenutation and a hour later it was still running

So, to summarise: what methods would you recommend for efficiently performing a signal split like this on uniformly sampled data? What if the signal were not uniformly sampled (i.e. a proper time-series) - fit and resample or...?

And why did my hf part still have the lf component?

Answer 1

The key question in the original post is: "why did my hf part still have the lf component?" since the answer to this dictates the rest.

The answer is that for a pure real signal, the Fourier coefficients are symmetric, i.e. if there are $n$ coefficients, the $k^{th}$ Fourier coefficient $F_k = F_{n-k}$. If the OP had taken a little longer to think about the theory before rushing to implement a (literally) half-baked idea, he might have seen this.

The original method to split a signal into lf and hf parts neglected the symmetry in Fourier coefficients for real signals and merely cleared the coefficients at one end before applying the inverse fourier & taking of real parts to recover the time domain signal.

rPad = PadRight[Take[fourierData,{1,splitF}],signalLen];
lPad = PadLeft[Take[fourierData,{splitF + 1,signalLen}],signalLen];
lfIF = InverseFourier[rPad,FourierParameters->{fourierA, fourierB}];
hfIF = InverseFourier[lPad,FourierParameters->{fourierA, fourierB}];
lfPart = Chop[Re[lfIF]];
hfPart = Chop[Re[hfIF]];

Thus in taking or clearing at only one end of the signal, only half the necessary coefficients were being used, hence the fact that in the illustrations provided, the spectrum of the lf part shows that the desired 261Hz signal has exactly half the amplitude of the 261Hz component in the input.

The correct approach is to replace the left/right padding with left-right/centre clearing

lfPad = ReplacePart[fourierData,Table[{i},{i,splitF+1,signalLen-splitF}]->0];
hfPad = fourierData -lfPad;
lfIF = InverseFourier[lfPad,FourierParameters->{fourierA, fourierB}];
hfIF = InverseFourier[hfPad,FourierParameters->{fourierA, fourierB}];
lfPart = Chop[Re[lfIF]];
hfPart = Chop[Re[hfIF]];

Note that given the symmetry on the Fourier coefficients, only the left hand parts of the spectra are plotted in the illustrations below (as was done by the OP, clearly demonstating that the OP was in fact aware of the symmetry but couldn't be bothered to think about it.)

Applying the revised code to a white noise signal shows that an effective brickwall filter is obtained, with the caveat that the signal can only be split at integral frequencies