Thanks everybody for the replies, but maybe I should explain a little
more!
The code I use to generate the test signal, is as follows;
void testSignal(buffer<complex<double>> *src)
{
double PI = 3.14159265359;
double incr = 1.0 / src->length;
double offset = 0.0;
for(long i=0L;i<src->length;i++)
{
double term1 = sin(40.0 * PI * offset) * exp(-100.0 * PI * pow(offset -
0.2,2.0));
double term2 = (sin(40.0 * PI * offset) + 2. * cos(160.0 * PI * offset))
* exp(-50.0 * PI * pow(offset-0.5,2.0));
double term3 = 2.0 * sin(160.0 * PI * offset) * exp(-100.0 * PI *
pow(offset - 0.8,2.0));
double t = (term1 + term2 + term3);
src->data[i].real(t);
offset += incr;
}
}
This basically generates three waveforms, the first and third are unique,
and the middle wave is a combination of the first and third (I swiped the
method for generating the signal from a wavelet book!)
my filters are based upon the daubechie (is that the correct spelling?)
4-tap co-efficients, which are computing using this function;
void reset()
{
double t = 4.0 * sqrt(2.0);
h[0] = (1.0 + sqrt(3.0)) / t;
h[1] = (3.0 + sqrt(3.0)) / t;
h[2] = (3.0 - sqrt(3.0)) / t;
h[3] = (1.0 - sqrt(3.0)) / t;
g[0] = h[3]; g[1] = -h[2]; g[2] = h[1]; g[3] = -h[0];
}
I generate the filter by implementing the recursion algorithm, on these
co-efficients.
When I perform the steps I described earlier, the final result, looks sort
of how I expected, being that its almost correct, apart from the fact that
the outputted waveform of the filtered test signal, in the time domain is
quite obviously symmetrical. i.e. the first half, and second half are a
mirror of themselves. I wasn't sure whether this was down to coding
error, or the results were suppose to look like this (as a result of the
chosen test signal), hence the need for another computer program to
compare the results with!
My background is software engineering rather than audio processing and I'm
mixing this work with all the other stuff I'm suppose to be doing at the
moment ;)
Regards Phill.
>Jerry Avins wrote:
>
>> Ben Bridgwater wrote:
>>
>> ...
>>
>>> http://ccrma.stanford.edu/~jos/parshl/Overlap_Add_Synthesis.html
>>> http://ccrma.stanford.edu/~jos/parshl/Choice_Hop_Size.html
>>>
>>> Is the theory here incorrect?
>>
>>
>> No, but you misinterpreted it. By "window", It seems that the author
>> means the desired filter function.
>
>Well, I think I did mistunderstand it, but in a different way... On
>re-reading it appears that he is referring to the window function, but
>that the system being desribed isn't actually filtering. It seems that
>he's using a windowed FFT with optimal overlap to smoothly sample the
>input, but then is simply using the input peaks to drive a sinewave
>synthesis output using overlap and add. It makes more sense now!
>
>Ben
>
FFT Based Filtering?
Started by ●November 20, 2006
Reply by ●November 22, 20062006-11-22
Reply by ●November 22, 20062006-11-22
On Tue, 21 Nov 2006 20:16:59 -0500, Ben Bridgwater <bbridgwater@myrealbox.com> wrote: (snipped)> >I've bookmarked that link, and intend to order a copy. > >Ben BridgwaterHi Ben, if you're a member of the IEEE Signal Processing Society, you can access Mark Borgerding's "filtering in the freq domain" article that he wrote for the March 2006 "DSP Tips & Tricks" column of the IEEE Sig. Proc. Magazine. His article's at: http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1598092 Good Luck, [-Rick-]
Reply by ●November 22, 20062006-11-22
moosedude wrote:> Thanks everybody for the replies, but maybe I should explain a little > more! > > The code I use to generate the test signal, is as follows; > > void testSignal(buffer<complex<double>> *src) > { > double PI = 3.14159265359; > double incr = 1.0 / src->length; > double offset = 0.0; > > for(long i=0L;i<src->length;i++) > { > double term1 = sin(40.0 * PI * offset) * exp(-100.0 * PI * pow(offset - > 0.2,2.0)); > double term2 = (sin(40.0 * PI * offset) + 2. * cos(160.0 * PI * offset)) > * exp(-50.0 * PI * pow(offset-0.5,2.0)); > double term3 = 2.0 * sin(160.0 * PI * offset) * exp(-100.0 * PI * > pow(offset - 0.8,2.0)); > > double t = (term1 + term2 + term3); > src->data[i].real(t); > > offset += incr; > } > } > > This basically generates three waveforms, the first and third are unique, > and the middle wave is a combination of the first and third (I swiped the > method for generating the signal from a wavelet book!) > > my filters are based upon the daubechie (is that the correct spelling?) > 4-tap co-efficients, which are computing using this function; > > void reset() > { > double t = 4.0 * sqrt(2.0); > > h[0] = (1.0 + sqrt(3.0)) / t; > h[1] = (3.0 + sqrt(3.0)) / t; > h[2] = (3.0 - sqrt(3.0)) / t; > h[3] = (1.0 - sqrt(3.0)) / t; > > g[0] = h[3]; g[1] = -h[2]; g[2] = h[1]; g[3] = -h[0]; > } > > I generate the filter by implementing the recursion algorithm, on these > co-efficients. > > When I perform the steps I described earlier, the final result, looks sort > of how I expected, being that its almost correct, apart from the fact that > the outputted waveform of the filtered test signal, in the time domain is > quite obviously symmetrical. i.e. the first half, and second half are a > mirror of themselves. I wasn't sure whether this was down to coding > error, or the results were suppose to look like this (as a result of the > chosen test signal), hence the need for another computer program to > compare the results with! >Hmmm. If the time domain is mirrored as you describe, that suggests that the frequency domain information applied to the IFFT is real. It should be complex with conjugate symmetry. John
Reply by ●November 22, 20062006-11-22
> >moosedude wrote: >> Thanks everybody for the replies, but maybe I should explain a little >> more! >> >> The code I use to generate the test signal, is as follows; >> >> void testSignal(buffer<complex<double>> *src) >> { >> double PI = 3.14159265359; >> double incr = 1.0 / src->length; >> double offset = 0.0; >> >> for(long i=0L;i<src->length;i++) >> { >> double term1 = sin(40.0 * PI * offset) * exp(-100.0 * PI * pow(offset->> 0.2,2.0)); >> double term2 = (sin(40.0 * PI * offset) + 2. * cos(160.0 * PI *offset))>> * exp(-50.0 * PI * pow(offset-0.5,2.0)); >> double term3 = 2.0 * sin(160.0 * PI * offset) * exp(-100.0 * PI * >> pow(offset - 0.8,2.0)); >> >> double t = (term1 + term2 + term3); >> src->data[i].real(t); >> >> offset += incr; >> } >> } >> >> This basically generates three waveforms, the first and third areunique,>> and the middle wave is a combination of the first and third (I swipedthe>> method for generating the signal from a wavelet book!) >> >> my filters are based upon the daubechie (is that the correctspelling?)>> 4-tap co-efficients, which are computing using this function; >> >> void reset() >> { >> double t = 4.0 * sqrt(2.0); >> >> h[0] = (1.0 + sqrt(3.0)) / t; >> h[1] = (3.0 + sqrt(3.0)) / t; >> h[2] = (3.0 - sqrt(3.0)) / t; >> h[3] = (1.0 - sqrt(3.0)) / t; >> >> g[0] = h[3]; g[1] = -h[2]; g[2] = h[1]; g[3] = -h[0]; >> } >> >> I generate the filter by implementing the recursion algorithm, onthese>> co-efficients. >> >> When I perform the steps I described earlier, the final result, lookssort>> of how I expected, being that its almost correct, apart from the factthat>> the outputted waveform of the filtered test signal, in the time domainis>> quite obviously symmetrical. i.e. the first half, and second half area>> mirror of themselves. I wasn't sure whether this was down to coding >> error, or the results were suppose to look like this (as a result ofthe>> chosen test signal), hence the need for another computer program to >> compare the results with! >> > >Hmmm. If the time domain is mirrored as you describe, that suggests >that the frequency domain information applied to the IFFT is real. It >should be complex with conjugate symmetry. > >John > >Hi John, can you expand on this a bit further? When I'm filtering in the frequency domain, I'm not quite sure what to do with the imaginary part of the data.. multiplying the imaginary part with the filter fft and the test signal fft, seemed to do nothing, as did leaving the imag. part alone, leaving it unprocessed. - Phill.
Reply by ●November 22, 20062006-11-22
moosedude wrote:> > > Hi John, > > can you expand on this a bit further? When I'm filtering in the frequency > domain, I'm not quite sure what to do with the imaginary part of the data.. > multiplying the imaginary part with the filter fft and the test signal fft, > seemed to do nothing, as did leaving the imag. part alone, leaving it > unprocessed. > > - Phill.Here's the idea. Matlab comes with a function like this, but I think this shorter m-file is easier to understand than their version. John ***************************** % fast convolution filter using overlap-add technique % FFT size is twice filter length function y = fastfilter(x,h) x = x(:); h = h(:); Nfft = 2*length(h); H = fft(h, Nfft); Npoints = Nfft - length(h); ip = 1; op = 1; y = zeros(2*length(x),1); while(ip+Npoints-1 <= length(x)) % grab an input block of length Npoints b = x(ip:(ip+Npoints-1)); % multiply FFT of input block by H, take inverse FFT B = fft(b,Nfft); B = B.*H; t = ifft(B); % concatenate all Nfft points of IFFT result onto output vector y(op:(op+Nfft-1)) = y(op:(op+Nfft-1)) + t; % advance input and output pointers by Npoints ip = ip + Npoints; op = op + Npoints; end y = y(1:(op+Nfft-1)); return
Reply by ●November 22, 20062006-11-22
"moosedude" <phill_holland@hotmail.com> writes:>> >>moosedude wrote: >>> Thanks everybody for the replies, but maybe I should explain a little >>> more! >>> >>> The code I use to generate the test signal, is as follows; >>> >>> void testSignal(buffer<complex<double>> *src) >>> { >>> double PI = 3.14159265359; >>> double incr = 1.0 / src->length; >>> double offset = 0.0; >>> >>> for(long i=0L;i<src->length;i++) >>> { >>> double term1 = sin(40.0 * PI * offset) * exp(-100.0 * PI * pow(offset > - >>> 0.2,2.0)); >>> double term2 = (sin(40.0 * PI * offset) + 2. * cos(160.0 * PI * > offset)) >>> * exp(-50.0 * PI * pow(offset-0.5,2.0)); >>> double term3 = 2.0 * sin(160.0 * PI * offset) * exp(-100.0 * PI * >>> pow(offset - 0.8,2.0)); >>> >>> double t = (term1 + term2 + term3); >>> src->data[i].real(t); >>> >>> offset += incr; >>> } >>> } >>> >>> This basically generates three waveforms, the first and third are > unique, >>> and the middle wave is a combination of the first and third (I swiped > the >>> method for generating the signal from a wavelet book!) >>> >>> my filters are based upon the daubechie (is that the correct > spelling?) >>> 4-tap co-efficients, which are computing using this function; >>> >>> void reset() >>> { >>> double t = 4.0 * sqrt(2.0); >>> >>> h[0] = (1.0 + sqrt(3.0)) / t; >>> h[1] = (3.0 + sqrt(3.0)) / t; >>> h[2] = (3.0 - sqrt(3.0)) / t; >>> h[3] = (1.0 - sqrt(3.0)) / t; >>> >>> g[0] = h[3]; g[1] = -h[2]; g[2] = h[1]; g[3] = -h[0]; >>> } >>> >>> I generate the filter by implementing the recursion algorithm, on > these >>> co-efficients. >>> >>> When I perform the steps I described earlier, the final result, looks > sort >>> of how I expected, being that its almost correct, apart from the fact > that >>> the outputted waveform of the filtered test signal, in the time domain > is >>> quite obviously symmetrical. i.e. the first half, and second half are > a >>> mirror of themselves. I wasn't sure whether this was down to coding >>> error, or the results were suppose to look like this (as a result of > the >>> chosen test signal), hence the need for another computer program to >>> compare the results with! >>> >> >>Hmmm. If the time domain is mirrored as you describe, that suggests >>that the frequency domain information applied to the IFFT is real. It >>should be complex with conjugate symmetry. >> >>John >> >> > > Hi John, > > can you expand on this a bit further? When I'm filtering in the frequency > domain, I'm not quite sure what to do with the imaginary part of the data.. > multiplying the imaginary part with the filter fft and the test signal fft, > seemed to do nothing, as did leaving the imag. part alone, leaving it > unprocessed.Hey Phill, You can see a homework report I did which performs the overlap/save method at http://www.digitalsignallabs.com/ece713/hw/hw4/hw.zip http://www.digitalsignallabs.com/ece713/hw/hw4/hw.pdf The .zip file includes the Matlab code. -- % Randy Yates % "How's life on earth? %% Fuquay-Varina, NC % ... What is it worth?" %%% 919-577-9882 % 'Mission (A World Record)', %%%% <yates@ieee.org> % *A New World Record*, ELO http://home.earthlink.net/~yatescr
Reply by ●November 22, 20062006-11-22
moosedude skrev:> Hello, > > I'm writing a simple C++ application, which performs the following tasks; > > 1) generates high/low pass daub 4-tap filters > > 2) performs a fft of those filters > > 3) performs a fft of the input, test signal > > 4) multiples the fft of the filters, with the fft of the test signal > > 5) performs a ifft of the result of the previous step, to output the > filtered result of the test signal.Here's a detail I didn't notice before: Items 4) and 5) combined are abit ambiguous. Exactly what do you do here? The phrasings above give the impression that you do something like (Matlab notation): %%%%%%%%%%%%%%%%%%%%%%%%%%%% X = fft(x); H1 = fft(h1,length(x)) ; % FFT of the first filter function H2 = fft(h2,length(x)); % FFT of the second filter function Y = X.*H1.*H2; % This is what item 4) says you do y = real(ifft(Y)); %%%%%%%%%%%%%%%%%%%%%%%%%%%% A better way would be to apply EITHER h1 OR h2 to filter the signal: %%%%%%%%%%%%%%%%%%%%%%%%%%%% X = fft(x); H1 = fft(h1,length(x)) ; % FFT of the first filter function H2 = fft(h2,length(x)); % FFT of the second filter function Y1 = X.*H1; % Apply LP filter Y2 = X.*Y2; & Apply HP filter y1 = real(ifft(Y1)); % LP-filtered data y2 = real(ifft(Y2)); % HP filtered data %%%%%%%%%%%%%%%%%%%%%%%%%%%% Rune
Reply by ●November 23, 20062006-11-23
Ben Bridgwater wrote:> Andor wrote:...> > No. FFT based filtering makes use of the fact that multiplication in > > frequency domain of two DFTed data frames (or vectors, as you prefer) > > corresponds to circular convolution in the time domain. All you need is > > zero-padding in the correct places to transform the circular into > > linear-convolution. No windows. > >OK, thanks. I had been thinking that the zero padding was only serving > the purpose of making sure that the entire impulse response of the > filter function was being included; I hadn't realized that it was also > serving to take away the circularity (which is why I thought the window > function was needed). Is there any intuitive way of seeing how this > works, or is it just a matter of falling out in the math?Well, the convolution y of two finite vectors, say h and x, with N1=length(h) and N2=length(x), is defined as y(n) = (x * h) (n) = sum_{k=0}^n x(k) h(n-k), and you can imediately see, that in general, y(n) =/= for n < N1 + N2, and y(n) == 0 for n >= N1 + N2. So if you want to filter in the frequency domain, the frame has to be long enough to accomodate y, which has length N1 + N2 - 1, hence the zero-padding. Regards, Andor
Reply by ●November 23, 20062006-11-23
> >moosedude skrev: >> Hello, >> >> I'm writing a simple C++ application, which performs the followingtasks;>> >> 1) generates high/low pass daub 4-tap filters >> >> 2) performs a fft of those filters >> >> 3) performs a fft of the input, test signal >> >> 4) multiples the fft of the filters, with the fft of the test signal >> >> 5) performs a ifft of the result of the previous step, to output the >> filtered result of the test signal. > >Here's a detail I didn't notice before: > >Items 4) and 5) combined are abit ambiguous. Exactly what >do you do here? The phrasings above give the impression >that you do something like (Matlab notation): > >%%%%%%%%%%%%%%%%%%%%%%%%%%%% >X = fft(x); >H1 = fft(h1,length(x)) ; % FFT of the first filter function >H2 = fft(h2,length(x)); % FFT of the second filter function > >Y = X.*H1.*H2; % This is what item 4) says you do > >y = real(ifft(Y)); >%%%%%%%%%%%%%%%%%%%%%%%%%%%% > >A better way would be to apply EITHER h1 OR h2 to filter >the signal: > >%%%%%%%%%%%%%%%%%%%%%%%%%%%% >X = fft(x); >H1 = fft(h1,length(x)) ; % FFT of the first filter function >H2 = fft(h2,length(x)); % FFT of the second filter function > >Y1 = X.*H1; % Apply LP filter >Y2 = X.*Y2; & Apply HP filter > >y1 = real(ifft(Y1)); % LP-filtered data >y2 = real(ifft(Y2)); % HP filtered data >%%%%%%%%%%%%%%%%%%%%%%%%%%%% > >Rune > >Hi Rune, you are correct, statement 4 & 5 are completely ambiguous, I am applying, as you suggested, the LP and HP separately on the same piece of data. Its part of the overall algorithm I'm applying to my source data, called a DWT (discrete wavelet transform). I may have asked some questions in the past here (a few months ago), relating to the same subject. The idea is that I'm an eventually going to adapt this DWT algorithm, to perform another algorithm called the CWT. (continuous wavelet transform) I must say, I'm overwhelmed by the amount of replies I've had to my initial question, I have never visited such a helpful forum before! I'm really appreciative to everyone for the response. I think, the problem is now solved, and as usual it was down to something stupid, when I was multiplying the fft results together, I was performing a standard arithmetic multiplication rather than a complex number multiplication !! oops. Regards Phill.
Reply by ●November 23, 20062006-11-23
moosedude wrote: ...> Its part of the overall algorithm I'm applying to my source data, called a > DWT (discrete wavelet transform). I may have asked some questions in the > past here (a few months ago), relating to the same subject. The idea is > that I'm an eventually going to adapt this DWT algorithm, to perform > another algorithm called the CWT. (continuous wavelet transform)Phil, I try neither to think much about transforms or to write about them here, but you just opened a new mystery for me. How can one do anything continuous on a finite computer. Your data are samples, are they not?> I must say, I'm overwhelmed by the amount of replies I've had to my > initial question, I have never visited such a helpful forum before! I'm > really appreciative to everyone for the response.Interesting problems often come from interesting people.> I think, the problem is now solved, and as usual it was down to something > stupid, when I was multiplying the fft results together, I was performing > a standard arithmetic multiplication rather than a complex number > multiplication !! oops.The official term for that is "Forrest-trees Syndrome". It happens to many explorers. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
