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?
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$