Hello everybody... I am trying to do the following but having troubles. The following is pseudo code.

(This happens once)
-Zero pad a 511 sized array(this will be used for overlap)

(The following lines loop)
-Capture 512 samples of an audio signal and store them in the first half of a 1024 size array
-Zero pad the remaining half of the array
-Window the first 512 samples with Hanning
-Convert the entire array to the frequency domain with FFT
-Multiply the contents with the frequency of an Impulse response
-Inverse FFT to get the filtered audio signal
-Add the overlap to the filtered signal
-Strore last 511 samples of the filtered signal to the overlap array
-Output the filtered signal
-Go back to capturing 512 samples etc...

The result is a filtered signal, with no clicks or pops (thanks to the windowing), but there are rapid amplitude chages. I have read just about every note I could find on the web, gone though PDF books and guides...but have not quite understood how many windowing do I do? Do I do another windowing operation somewhere? There are mentions online in places about how its done, but some just suggest it too briefly, and some go into some deep DSP theories which I admit is not one of my strong points. Could somebody who has an understanding on this topic help me out? Pseudo code or c/C++ outlines would be great.

Thanks a lot in advance!

PS: This is for a suit of plug-ins I am writing, whcih will be freely available at the end of the month. I will let you know more about them and where you can download them. But first I need to make sure they work

I think you should try this approach :

Construct an array of N samples with N/2 samples shifted from previous frame.. DO NOT zero pad any samples !!

Windowed the array with a suitable window..

Pass the windowed array into a FFT


At the reconstruction side, the IFFT would reproduce the same windowed array.. Overlapped and add with the previous half window..

eg :
w(n)S_2(n) + w(N/2 + n)S_1(N/2 +n) = S_Original(n)..

Since that S_2(n) = S_1(N/2 +n) = S_Original(n) ;

then the windowing function MUST obey this equation for perfect reconstruction.. :

w(n) + w(N/2 +n) = 1


One possible window solution is the hanning (or was it hanming ???) window given by the following equation :

w(n) = 0.5*(1.0 - cos((2*PI * (n+0.5) / N))) /// Sorry for the mistake !!!

It seems like you want to do plain filtering by "fast convolution" (via FFT).
There's NO need to apply a window function !

Let 'k' be the length of the impulse response of the filter.
Now, if you want an infinite signal to be convolved with this impulse response you can do this block-wise and do some overlap adding. This is legal because of the properties of the filter (being a time-invariant linear system)

Let 'x' be your signal you want to be convolved with the filter's IR.
Let 'n' be the smallest power of 2 which is graeter than or equal to k-1.

You can take blocks of n samples of your signal, convolve them seperately with the impulse response and do the overlap-add of the resulting signal blocks which have at maximum n+k-1 non-zero-samples (due to the convolving).
because of k-1<=n we know n+k-1<=2n. So, the length in samples of the convolved signal blocks are 2n at maximum. This convolution can be done via the FFT of size 2n.

Here's an example:


The zero-padding is actually needed because the convolution you can calculate via FFT is a cyclic one and our source signal is NOT periodic.

HTH,
Sebastian

(1) It's called "Hann" window. See Hann and Hamming Window

(2) Your approach is flawed and introduces temporal alias artefacts. It does not take into account that FFT convolution is cyclic. (You need zero padding to compensate for this)

Sebastian

Thank you for your replies and suggestions. I got rid of any windowing and taken Sebastian's sugegstion. I have been advised by others to do just that and aparrently it works! I have adjusted my code to the sugegsted method (which actually is what it was before I got carried away with the Windowing theory). However I get clicks/pops at the end of each frame block processed. This results in buzzing. I have added the code of my processing with comments bellow in case anybody can spot a problem.

Regarding the spectrum multiplication. I am using Apple's vDSP library's function zvmulD(). In the documentation it is described as "Complex Vector Multiply." On of the function arguments is the "conjugate." The documentation says "Assign conjugate flag 'conjugate' a value of 1 for normal multiplication or -1 for multiplication by conjugated values of input 1." I have tried both, and both give different waveforms and both still have the clicks/pops. I am new to FFTs and DSP in general so is it normal multiplication I need for my scenario or multiplication by conjugated values of input 1? While both give clicks and pops, perhaps this may shed some light on the overall problem with the pops in the signal.


Thanks a lot!

YES! Its solved! The problem was the size of the array that the complex vector multiplication was being done for. It was log2n which in this case equals 10, which it should have been nOver2 which equals 512. I am putting the correct code bellow. It should be useful to newbies like myself who learn DSP better by reading code than text and math formulas.

Thanks for your help!!!

I guess it's correct.

you may want to change the last part to

which saves some assignments.
Make sure that the overLap buffer contains 512 samples.

BTW: With a complex-valued FFT you can convolve a stereo signal (treat is as quadrature signal) with a mono (real-valued) impulse response at once for both channels or convolve a mono (real-valued) signal with a stereo (quadrature) impulse-response for both ears in case of the HRTF stuff. There's no need to do this processing you are doing twice for each channel.

bye,
Sebastian

Really? What is this temporal alias artefacts ?

But isn't that padding the 2nd half of the FFT input samples with zero is equivalent to the multiplication with a window of a "brickwall" shape / the rectangular window ?

As mentioned.. the window must obey the equation : w(N/2+n) + w(n) = 1 which is also true here.. But this window shape could produce boundary artifacts for the reconstructed audio as it is obvious in cases of DCT used in image and video compression systems..

Windowing the input samples with an "appropriate window function" would spread the boundary blocking artifacts into 2 adjacent frames.. and kept well below the masking properties of the actual signal itself..


wkwai

If you convolve a discrete signal which consists of n consecutive non-zero samples with another signal with k samples you get a signal with n+k-1 consecutive non-zero samples (or less). If you try to compute this via an n-point FFT the last k-1 samples of the convolved signal which should contribute to the next block get wrapped around and create an "alias". That's why I used the term "cyclic convolution". If you want to calculate this via an n-point FFT you need to shorten the blocks you convolve with the impulse response so the convolved blocks are at most n samples long. Your idea only works in case of the impulse response being a scaled unit impulse where k=1.

recap:
b = blocksize / step
n = FFT-Size
k = length of the impulse response in samples
Then this must be satisfied:
b+k-1<=n
The signal blocks and the impulse response have to be zero-padded to n samples.

FFT-Convolution-Example:
In case of n=4 and the two signals to be convolved via the FFT are a=(a_0,a_1,a_2,a_3) and b=(b_0,b_1,b_2,b_3) the "FFT-Convolution" yields c=(c_0,c_1,c_2,c_3) with
c_0 = a_0 * b_0 + a_3 * b_1 + a_2 * b_2 + a_1 * b_3
c_0 = a_1 * b_0 + a_0 * b_1 + a_3 * b_2 + a_2 * b_3
c_0 = a_2 * b_0 + a_1 * b_1 + a_0 * b_2 + a_3 * b_3
c_0 = a_3 * b_0 + a_2 * b_1 + a_1 * b_2 + a_0 * b_3

Notice when b_1>0, a_3 contributes to c_0 which should contribute to c_4 instead in case of a "normal" convolution. This is a cyclic convolution.


HTH,
Sebastian

Remember that a convolution is linear and time-invariant. It's therefore legal to split the signal into blocks, convolve them seperately and add these convolved blocks to form the new signal.

Try to understand what actually happens if you multiply 2 Fourier spectra and why I always used the term cyclic convolution. You may also want to check Chapter 18 of "DSPGUIDE".

Sebastian

Are you suggesting that if I input an N samples into a FFT and then the output of the FFT back into a IFFT, I don't get back the original N samples ??? assuming that the normalization factors are properly handled and truncation errors kept to minimum..

I will try to write a simple program to verify this.. Thanks for the "enlightenment"..

wkwai

No.

Let a,b \in \mathbb{C}^n

where ".*" denotes the component-wise product. To avoid those nasty wrap-arounds 'a' and 'b' have to be zero-padded. Be 't' the greatest index with a_t!=0 and 'k' the greatest index so that b_k!=0. If k+t-1<=n there won't be any wrap-arounds.

I don't know how to explain this any clearer. Keep in mind that the Fourier transform gives you a representation of a periodic signal.

For further studies I suggest reading some DSPGUIDE chapters.

Sebastian

I believe the simple interpretation of "linear convolution in the time domain = multiplication in the frequency domain" only applies to continuous time/frequency signals, where you don't have sampling occuring and things wrapping around because of it. For the case of the FFT, which is a transform based on discrete time, the equivalent of multiplication in frequency is circular/cyclic convolution in the time domain, rather than your typical linear convolution. ie. FFT-based convolution is the same as convolving your discrete time signal wrapped around itself, which gives you a periodic signal. Thus you need to zero pad it enough to increase the period of the signal in avoid the effects of circular convolution ie. the effect of terms that have wrapped around

In subband image coding, we usually don't use zero padding as the discontinuity can cause edge distortions, but opt for something smoother like symmetric extension instead.

Here's a matlab example:



Now that's linear convolution. Doing a straight multiplication using FFT gives:



which isnt the same.

Now we do some zero padding. Since the final output is 6+3-1=8 samples, I zero pad them out to 8:



Now we have our linear convolution result

Yeah.. I think I understand the problem here.. We are dealing with short time convolution here.. Not continous time convolution.. As a result, the extra k-1 samples need to be regenerated from the IMDCT..

I did some studying on the short time convolution.. and verified that the k-1 extra samples belong to the next frame.. which you need to add to first few samples of the next frame..

In my ealier equations, I assumed that convolution is just a straight forward multiplication in the freq domain.. I failed to take into consideration of the "boundary conditions" needed to perfectly reconstruct the audio signal..

My appologies..

But I have another very troubling problem.. Supposed that I had an FFT spectral of N/2 spectrals components.. Then I attempt to zeroed out the upper coefficient to "band-limit" the spectral.. According to this insight provided from you folks.. this is WRONG!!!

What about the CoolEdit FFT filter option where one could just generate "ANY" filter frequency response shape as he likes?? I wondered how is it implemented?
Did they actually reconvert this response drawn by any user into impulse response.. by an IFFT and then determined the input samples length to the FFT ??

I wondered too if the actual CoolEdit FFT filter response is just an approximation of the one drawn by the user??

In Twin VQ / VQF, there is also some very troubling problem.. because, Mr. Arigatou of NTT who developed the coding technique assumed that the MDCT spectral shape is identical to the FFT spectral shape.. So, the estimated LPC spectral shape is used to flattened the MDCT spectral shape which is not correct.. mathematically.. since the lpc calculation is calculated in the TIME DOMAIN and the MDCT decomposition is actually a critically downsampled subband filterbank !!

Why shouldn't the convolution be done in time domain instead of trying to transform everything into freq domain.. and then having to calculate the square 2 root of the LPC estimated spectral shape ?? It is an unnecessary algorithmic complexity.. plus the fact that it is mathematically incorrect..



wkwai

I don't know for sure, but I guess Cool Edit generates an impulse response of the desired filter's frequency response via iFFT + windowing and uses it for fast convolution via FFT the way I described above.


Yes, but a close one. It depends on your desired frequency response and FFT size settings. If you want a very close approximation of a brickwall lowpass you need the impulse response to be vera large.


Well, the thing is: You don't convert it back to time-domain after frequency domain filtering, so it's totally reversible. The actual signal itself isn't altered after encoding and decoding. But the introduced noise get's "filtered" through the LPC synthesis filter which is done in the MDCT domain and introduces an temporal alias. But it does not matter because it's only noise.

Sebastian

[/quote]
Well, the thing is: You don't convert it back to time-domain after frequency domain filtering, so it's totally reversible. The actual signal itself isn't altered after encoding and decoding. But the introduced noise get's "filtered" through the LPC synthesis filter which is done in the MDCT domain and introduces an temporal alias. But it does not matter because it's only noise.

Sebastian

[/quote]


Sure it is not a problem at the encoder-decoder pair.. but at the encoder itself.. It is assuming that the mdct spectral shape is FFT spectral shape !!! Which is mathematically wrong.. Would have been better to convolute in time domain first before transforming it with the MDCT.. ??? That would make something like a GSM speech lpc analysis.. time-domain..

wkwai

The MDCT is related in the same way as the FFT to the LPC analysis/synthesis.

You could have done it in the time domain. But in this case it would have been very expensive because of the high orders of the filters. In fact you need higher orders if you want to do this in the time-domain since you cannot use the logarithmic / bark frequency scale like it's done in VQF / Vorbis (with floor 0).


Sebastian

I also suspect that is what the CoolEdit FFT filter option does. It can use either the window method or frequency-sampling method of FIR filter design.

Perhaps.. it is in their (Japanese) Engineering mindset.. Engineering is both an art form as well as Science ?? I know that in some Engineering cultures.. Engineering students are taught to design like an "artist".. which is unfortunately.. I am NOT trained in.. which is also why I DON'T do things like GUI designs..