0
$\begingroup$

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.

enter image description here

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?

$\endgroup$
  • $\begingroup$ Why not use LowpassFilter and HighpassFilter? There is no such thing as an "ideal" filter, and with these commands you can get closer by choosing the filter length option larger. $\endgroup$ – bill s Nov 23 '17 at 18:35
  • $\begingroup$ You may find some details on Fourier here. Have you taken into account the fact that the spectrum has the frequency range you expect (up to sample rate/2) followed by the values corresponding to negative frequencies? $\endgroup$ – Hugh Dec 3 '17 at 19:56
  • $\begingroup$ @Hugh That is indeed exactly the issue I identified in my answer below :) $\endgroup$ – Julian Moore Dec 4 '17 at 8:19
2
$\begingroup$

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.)

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ (+1) for the use of ReplacePart, i'll try to make this function part of my day-to-day repertoire now. $\endgroup$ – Thies Heidecke Dec 3 '17 at 16:02
  • $\begingroup$ @ThiesHeidecke That's nice - it's the first time any MMA function usage I have employed has been noteworthy to anyone! $\endgroup$ – Julian Moore Dec 4 '17 at 8:21

Your Answer

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

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