Normalization of zero padded signals

I have a simple question regarding zero padding and normalization. Consider an impulse resonse of a 4 point moving average filter. and its fft zero padded to 1024 points..

x=[1/4 1/4 1/4 1/4]
X=fft(x,1024 )
xpowrsum=dot(x,x)
Xpowrsum=dot(abs(X),abs(X))/1024
plot(fftshift(abs(X)))

By Parsevals theorem the two energies are equal as expected. However, the fft without scaling shows the correct frequency response with a gain of 1 at 0 Hz. So why do I always read the FFT should be scaled by the number of samples before zero padding (in this case 4) if I am interested in the magnitude response of the filter?

Answer 1

[Matt J]

So why do I always read the FFT should be scaled by the number of samples before zero padding (in this case 4) if I am interested in the magnitude response of the filter?

The FFT is a tool with many applications, each with its own appropriate scaling.

Scaling by 1/N is done when the FFT is being used to evaluate the Discrete Fourier Series.

When it is being used to approximate the continuous Fourier transform, it is scaled by the time sampling interval 1/Fs.

To achieve Parseval's equality, the fft should be scaled by 1/sqrt(N):

x=[1/4 1/4 1/4 1/4];
X=fft(x,1024 )/sqrt(1024);
xpowrsum=norm(x).^2

xpowrsum = 0.2500

Xpowrsum=norm(X).^2

Xpowrsum = 0.2500

[Marc Fuller]

Your computation of Parseval's theorem and the one in the original script are computationally equivalent. I am not approximating a continous signal. I I am trying to calculate the frequency response of a 4 point discrete time moving average filter. If you scale the resulting FFT by anything besides 1 the magnitude at 0 will no longer be one, which would be incorrect.

[Matt J]

Yes, it is equivalent, but you asked "why do I always read the FFT should be scaled by the number of samples ..." and I've explained to you that it shouldn't, unless you are trying to compute the Discrete Fourier Series. In a Discrete Fourier Series, the DC harmonic is the mean of the signal values, not the sum of the signal values.

If you think your DC value should be 1, for whatever it is that you're doing, then you are correctly scaling the FFT for whatever it is that you're doing.

[Paul]

@Matt J

Is it standard practice to include the 1/N factor in the definition of the DFS coefficients? I feel like I've also seen it the other way, i.e., the DFS coefficients don't include the 1/N, as is the case in fft, and the 1/N is then applied to the sum that reconsructs the signal from those DFS coefficients.

[Marc Fuller]

I can tell you that there is not a consistent definition. If you have three different text books, chances are good you will see three definitions. However if you don't scale the fft, then the energy calculation (Parseval's Theorem), no matter how you do it becomes a problem. Note that if you just do ifft(fft(x)); in Matlab, you do get the original function back, its just the scale in the frequency domain does not reflect the power in the original signal.

[Matt J]

The definition of the DFS will indeed vary from textbook to textbook. The bottom line is if you want your DC component to be the average of the signal values, you would divide fft() by N. Otherwise, DC will be the sum of the signal values.

[Matt J]

One example to motivate the 1/N factor is to consider a periodic signal like,

If the goal is to recover the coefficients of the sinusoidal terms (5 and 3), we can see in the following code that the 1/N is necessary.

N=10; 
n=(0:9)';
x=5+3*exp(1j*2*pi*n/N);

c=fft(x)/N

c = 
   5.0000 + 0.0000i
   3.0000 - 0.0000i
   0.0000 + 0.0000i
   0.0000 - 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
  -0.0000 + 0.0000i
   0.0000 + 0.0000i
  -0.0000 + 0.0000i
   0.0000 + 0.0000i

Answer 2

[Marc Fuller]

So after thinking about this, here is the reason you should not normalize a filter.

Consider y(t)=x(t) convolved with h(t)- where h(t) is some filter such as the moving average filter in this question.

If we want to be consistent, suppose we transform y(t) with normaiization to Y(f) by dividing by the length of the fft, N

Now if we transfrom x(t) convolved with h(t) We have X(f)H(f) by the convolution theorem. We only want to normalize X(f) to have the product identical with Y(f). since Y(f)/N= (X(f)H(f)/N.

[Matt J]

No, that doesn't explain it. You could just add a factor of 1/N to the convolution threorem, just as you did with parseval's theorem.

FFT(conv(x,h))=FFT(x)*FFT(h)/N

[Marc Fuller]

So my assumption is that you are treating all signals that you transform into the frequency domain identically. Yes, you could transform the filter, but then you would not be treating all signals you transform identically, since you would then not normalize X(f). I am trying to differentiate between signals and filters.

[Matt J]

I'm not suggesting here that you transform the filter. I'm suggestign that the convolution theorem include an additional scale factor of N. Errata: above I should instead have had,

FFT(conv(x,h))=N( FFT(x).*FFT(h))

Now you are free to normalize all the FFTs by 1/N and th equation will still be true.

[Paul]

The governing equation for the linear system is

DTFT(h[n]*x[n]) = H(f)X(f) where H(f) = DTFT(h[n]) and X(f) = DTFT(x[n]).

To approximate the RHS of this equation for finite duration sequences h and x using the DFT as implemented by fft, we'd have

H(f)X(f) evaluated at samples of f is the product DFT(h[n],N) .* DFT(x[n],N) where N >> max(numel(x),numel(h)).

If the DFT were implemented with the 1/N factor, we'd have to multiply the DFT product by N^2 to recover the correct magnitude of the DTFT.

[Matt J]

The FFT as currently implemented in Matlab satisfies,

FFT(conv(x,h))= FFT(x).*FFT(h)

If we divide by N^2 on both sides this becomes,

FFT(conv(x,h)/N)/N= FFT(x)/N.*FFT(h)/N

which means the equation is still true if the FFT() operator is replaced by FFT()/N and the conv() operator is replaced by conv()/N.

[Paul]

Kind of getting off the topic of normalization, but I don't think this equation can be true

FFT(conv(x,h))= FFT(x).*FFT(h)

numel(conv(x,h)) = numel(x) + numel(h) - 1

but the numel of the RHS is equal to numel(x) (and numel(h) since numel(x) == numel(h) for the .* to work).

Illustrating my comment above with some data.

rng(100);
h = rand(1,10);  % FIR
x = rand(1,15);  % input
z = conv(h,x);   % output
[Z,w] = freqz(z,1,linspace(0,2*pi,1024)); % DTFT of output
H = fft(h,128);  % zero-padded DFT of impulse response
X = fft(x,128);  % zero-padded DFT of input
ww = (0:127)/128*2*pi; 
figure  % plot gain and phase.
hold on
stem(ww,abs(H.*X))
plot(w,abs(Z))

figure
hold on
stem(ww,angle(H.*X))
plot(w,angle(Z))

[Matt J]

Kind of getting off the topic of normalization, but I don't think this equation can be true FFT(conv(x,h))= FFT(x).*FFT(h)

Here, I was somewhat abusing the notation "conv()" which is meant to mean cyclic convolution. I demonstrate the equivalence below using my own cyclic convolution routine which I've attached.

x=rand(1,5); h=rand(1,5);

fft(cyconv(x,h))

ans = 
   6.3779 + 0.0000i  -0.5432 - 0.1599i   0.3127 + 0.3201i   0.3127 - 0.3201i  -0.5432 + 0.1599i

fft(x).*fft(h)

ans = 
   6.3779 + 0.0000i  -0.5432 - 0.1599i   0.3127 + 0.3201i   0.3127 - 0.3201i  -0.5432 + 0.1599i

[Paul]

I thought that you probably meant that. I haven't looked at cyconv. Is it preferred over Matlab's cconv for some reason?

rng(100);
x=rand(1,5); h=rand(1,5);
fft(cconv(x,h,5))

ans = 
   4.8831 + 0.0000i   0.1012 + 0.2000i  -0.0803 + 0.7535i  -0.0803 - 0.7535i   0.1012 - 0.2000i

fft(x).*fft(h)

ans = 
   4.8831 + 0.0000i   0.1012 + 0.2000i  -0.0803 + 0.7535i  -0.0803 - 0.7535i   0.1012 - 0.2000i

[Matt J]

It requires no toolboxes and works on arrays of any dimension.