FFT question

I tried to run this experiment in MATLAB:


numberOfSamples = 80;
sampleBlock = randn(1,numberOfSamples);
counter = 1;
for numberOfZeros = numberOfSamples:2000
x = [sampleBlock zeros(1,numberOfZeros-numberOfSamples)] ;
y = filter(H,1,x);
G1 = abs(fft(y,2*length(x)));
G2 = abs(fft(x,2*length(x))).*abs(fft(H,2*length(x)));

numberOfZerosVector(counter) = numberOfZeros;
varianceOfError(counter) = var(G1-G2);

counter = counter + 1;
end
figure
plot(numberOfZerosVector,varianceOfError)

It looks like I have to zero pad the sampleBlock with
around 2000 zeros before G2 is approximately equal to
G1 ? (H has 4001 filter coefficients).

My general question is:

If we have two known time series x and y and we know that y is
obtained by filtering x with
some unknown FIR filter coefficients H of known length N, how
do we best approximate the spectrum of y as a product of the spectra
of x and H ?

On Jan 29, 5:25&#4294967295;pm, Mauritz Jameson <mjames2...@gmail.com> wrote:
> I tried to run this experiment in MATLAB:
>
> numberOfSamples = 80;
> sampleBlock = randn(1,numberOfSamples);
> counter = 1;
> for numberOfZeros = numberOfSamples:2000
> x = [sampleBlock zeros(1,numberOfZeros-numberOfSamples)] ;
> y = filter(H,1,x);
> G1 = abs(fft(y,2*length(x)));
> G2 = abs(fft(x,2*length(x))).*abs(fft(H,2*length(x)));
>
> numberOfZerosVector(counter) = numberOfZeros;
> varianceOfError(counter) = var(G1-G2);
>
> counter = counter + 1;
> end
> figure
> plot(numberOfZerosVector,varianceOfError)
>
> It looks like I have to zero pad the sampleBlock with
> around 2000 zeros before G2 is approximately equal to
> G1 ? (H has 4001 filter coefficients).
>
> My general question is:
>
> If we have two known time series x and y and we know that y is
> obtained by filtering x with
> some unknown FIR filter coefficients H of known length N, how
> do we best approximate the spectrum of y as a product of the spectra
> of x and H ?
The answer is in the question is it not?
Y(w)=X(w).H(w)

That's one answer. I was wondering if there is a better way of
approximating Y
when your FFT signal buffer is small but your unknown FIR filter might
have a lot of taps?

On Sun, 29 Jan 2012 13:37:27 -0800, Mauritz Jameson wrote:

> That's one answer. I was wondering if there is a better way of
> approximating Y
> when your FFT signal buffer is small but your unknown FIR filter might
> have a lot of taps?

You are playing fast and loose with terminology, probably because you are 
confused about the relationship between the discrete Fourier transform 
(the FFT is a way to _realize_ the DFT), and the discrete-_time_ Fourier 
transform.

The discrete Fourier transform either works on sampled data vectors of 
finite length, or on sampled data vectors that are one cycle of an 
infinitely long, perfectly repeating sequence (which one it is depends on 
your background, and is the subject of an infinitely long flame war on 
this group).  The discrete-time Fourier transform works on signals that 
are sampled in time and are infinitely long.

Each of the two transforms has a spectra: the DFT (FFT) has a spectrum 
that is of finite length and that comes in discrete samples.  The 
discrete time Fourier transform has a spectra that is finite length (2 pi 
radians), and is continuous (not discrete) over that length.

When you invoke Matlab's "filter" function, it is assuming that you are 
asking it to filter a finite length sample of an infinite length signal, 
and it truncates the result in time.  Taking the FFT of the resulting 
finite-length sample is an exact process (not an approximation), but the 
result is an approximation of a discrete-time Fourier transform, because 
the finite-length sample is but an approximation of the infinite-length 
signal that the 'filter' function (and your spectra) assume.

So, you have three choices:  Keep everything in the sampled-time, finite 
length domain, keep everything in the sampled-time, infinite length 
domain, or accept that you are making approximations when you cross from 
one domain to another, and take steps to understand the approximation 
you're making and to minimize the deleterious effects of doing so.

To keep everything in the sampled-time, finite length domain, you must 
use a filter function that performs a circular convolution on a finite-
length sample.  This will make it 'fit' with the FFT, but it won't mean 
much in real life (unless you're wondering what your filter will do to a 
periodic signal with the period you've given it).

You're really operating in the discrete-time domain with signals of 
infinite (or practically infinite) extent.  So to do the analysis 
"properly" you should use the discrete-time Fourier transform (which 
pretty much means that you need to solve everything symbolically).

What you're _trying_ to do is to take the third path, and use the FFT for 
a signal in the discrete-time (but not finite-extent) world.  
Fortunately, this is fairly easy to do if your signal and your filter are 
both of finite extent: take your signal and pad it with as many zeros as 
your FIR filter has taps, and take your FIR filter impulse response and 
pad it with as many zeros as your signal has duration (I'd pad both with 
extra zeros, to reassure myself that I didn't screw anything up, and 
maybe to bring things up to some 2^N length).  Now do your whole FFT 
thing the way you're doing it.  Because it's a _finite_ impulse response 
filter and a _finite_ signal, you can rest assured that the filter output 
will be finite, and no longer than the signal duration plus the filter 
length: any extra zero-padding would just get you more zeros.

Doing it this way, with adequate zero padding, will give you a 
theoretically exact answer: any mismatch between a convolution of filter 
and signal, vs. the ifft of the product of the spectra, will be due to 
finite numerical precision in the calculation, not fundamental 
theoretical errors in what you are doing.

Of course, if either your filter or your signal is of infinite extent, 
then you need to approximate them by truncation, and you're back to 
needing to understand what you're doing.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com