5
$\begingroup$

I have computed Fourier transform from ecg data.It does obey Parseval's theorem, relation

$$\sum_{n=0}^{N-1} \Big| x[n] \Big|^2 = \frac{1}{N} \sum_{k=0}^{N-1} \Big| X[k] \Big|^2$$

is fulfiled. But after I use window function (Blackman,Hann,Hanning) energy is not conserved anymore.Is this normal? If not,how can I normalize my windowed fft?

I calculate fft this way:

    N = signal.shape[0]
    frq = np.fft.fftfreq(N,freq)
    frq = frq[range(N/2)]
    window = window_func(N)
    fft_y = np.fft.fft(signal*window)/np.sqrt(N)
    fft_y = Y[range(N/2)]
    plt.plot(frq, abs(fft_y))

And thats how it look when I plot. I'm a little worried about how almost all windows drastically reduce amplitude,again is this normal? enter image description here

Thank you in advance for all help.

$\endgroup$
3
$\begingroup$

windowing your data, $x[n]$, with window $w[n]$ is equivalent to windowing the square of your data (or square-magnitude), $\Big|x[n]\Big|^2$, with the square of the window $w^2[n]$. so think of this square of the window $w^2[n]$ as just another window.

now stochastically your typical squared sample $x^2[n]$ will have some expected mean-square. the notation we use is either the left-hand or right-hand side:

$$ \overline{\Big|x[n]\Big|^2} = \mathscr{E}\Big\{ \Big|x[n]\Big|^2 \Big\} $$

if $x[n]$ was truly random data with a mean and variance, we can say that the mean-square is $ \overline{\Big|x[n]\Big|^2} $.

now this is only true if $x[n]$ is sufficiently random, in comparison to and independent of window $w[n]$:

$$\begin{align} \overline{\Big|w[n] \cdot x[n]\Big|^2} &= \overline{ w^2[n] \cdot \Big|x[n]\Big|^2 } \\ &= \overline{w^2[n]} \cdot \overline{ \Big|x[n]\Big|^2 } \\ \end{align}$$

then:

$$ \begin{align} \sum_{n=0}^{N-1} \overline{\Big|w[n] \cdot x[n]\Big|^2} &= \overline{w^2[n]} \cdot \sum_{n=0}^{N-1}\overline{ \Big|x[n]\Big|^2 } \\ &= \overline{w^2[n]} \cdot \frac{1}{N} \sum_{k=0}^{N-1} \Big| X[k] \Big|^2 \\ \end{align} $$

where

$$ \overline{w^2[n]} = \frac{1}{N} \sum_{n=0}^{N-1} w^2[n] $$

this scaling feature is based on the assumption of $x[n]$ being random and independent of the window $w[n]$. it's really then an approximation for your specific application that gets more accurate as the width of the window (which is $N$) gets large in comparison with any noticeable period of $x[n]$.

$\endgroup$
  • 1
    $\begingroup$ The scaling (normalization) coefficient can be computed accurately. There's no need for stochastic anything here :). I know that there was an article describing this in detail so I'll try to find it. $\endgroup$ – dsp_user Mar 6 '18 at 12:13
  • 1
    $\begingroup$ well, already i know that you're mistaken, mathematically. consider a goofy window $w[n]$ that happens to be zero when $x[n]$ is not zero, but $w[n]$ is non-zero for at least some times that $x[n]$ is zero. what's going the be the power of the product, $w[n]\cdot x[n]$? what will be the scaling coefficient you will derive from the window $w[n]$? zero ?? $\endgroup$ – robert bristow-johnson Mar 6 '18 at 19:28
  • $\begingroup$ I appreciate your comment. I never mentioned that what I wrote in my post is mathematically absolutely correct. I did use the word accurately in a comment and perhaps I shouldn't have. The article Jojek linked to is a way to compute normalized magnitude values and this is often useful. Many applications use similar approaches to calculate magnitudes but it seems that I failed to convey this adequately. In conclusion, if the quoted article is of any relevance (and there are people who think it is), we shouldn't dismiss it easily. $\endgroup$ – dsp_user Mar 6 '18 at 20:11
  • $\begingroup$ @dsp_user The point is that Jojek is answering a different question than the one posed here. $\endgroup$ – Jazzmaniac Mar 6 '18 at 20:15
8
$\begingroup$

Parseval's theorem will hold, but take into account that your signal in the time domain will no longer be $x[n]$. Namely, if you have that

$$\sum_{n=0}^{N-1} \Big| x[n] \Big|^2 = \frac{1}{N} \sum_{k=0}^{N-1} \Big| X[k] \Big|^2$$

then, if you window the signal $x[n]$ with a window $w[n]$, your signal will now be $\hat{x}[n]=x[n]w[n]$, and the theorem will now return the following equality:

$$\sum_{n=0}^{N-1} \Big| \hat{x}[n] \Big|^2 = \frac{1}{N} \sum_{k=0}^{N-1} \Big| \hat{X}[k] \Big|^2$$

where $\hat{X}[k]$ is the DFT of $\hat{x}[n]$ (and obviously differs from $X[k]$).

$\endgroup$
  • 1
    $\begingroup$ Yes, but it's still possible to account for the energy loss due to the windowing. The easiest way is to compute a windowing function coefficient /weight and then divide magnitudes by this coefficient. In this way, our magnitudes will reflect the original energies, which is often desirable (e.g a unit amplitude signal will correspond to 0 dB etc). $\endgroup$ – dsp_user Mar 6 '18 at 8:44
  • $\begingroup$ @dsp_user Hi, could you elaborate how to compute windowing function coefficient? I'm assume that for every window type there is its own coefficient, right? $\endgroup$ – wiedzminYo Mar 6 '18 at 10:01
  • $\begingroup$ weidzmin, I did that. just read the answer. $\endgroup$ – robert bristow-johnson Mar 6 '18 at 11:20
3
$\begingroup$

While the Fourier transform, discrete or continuous, can be regarded as unitary transform i.e a naturally norm preserving change between orthonormal bases in a normed complex vector space, the windowed FT does not in general possess this quality. And non-unitary operators cannot be turned into unitary ones by re-scaling. Here is why:

Unitary operators can be diagonalised over the complex numbers and all eigenvalues are of magnitude 1. That means if you put a vector of a certain norm in, you get a vector of that same norm out, no matter what the vector is. The Fourier transform operator $F$ is unitary.

Applying a window function can also be regarded as applying a linear operator $W$ to a signal. This operator is diagonal in the domain where the point wise multiplication of the window happened, so usually the time domain. The eigenvalues are then the individual values of the window function. These are usually not all the same unless you have a constant magnitude window. So in general the eigenvalues of the window operator cannot be rescaled to be of magnitude 1.

The composite linear operator of first applying the window function and then applying the Fourier transform cannot be easily diagonalised and the eigenvalues therefore not determined right away. However, we can say something about the norm of the vectors that go out and those that go in allowing us to make a statement about the eigenvalues of the composite operator.

If we take a window eigenvector $v_1$ with eigenvalue $\lambda_1$ as input, we know that the norm of the output is $\left|FW\,v_1\right|=|\lambda_1| \cdot |v_1|$. Applying this to the largest and smallest (by magnitude) eigenvalues, $\lambda_{\mathrm{min}}$ and $\lambda_{\mathrm{max}}$ respectively, of the window operator, we get a spectral range of $FW$ that must contain the subset $\{|\lambda_{\mathrm{min}}|,|\lambda_{\mathrm{max}}|\}$. The spectral range of $FW$ is bounded by the magnitudes of its own eigenvalues, so that we know that smallest and greatest eigenvalues (again, by magnitude), $\eta_{\mathrm{min}}$ and $\eta_{\mathrm{max}}$ respectively, are bounded: $|\eta_{\mathrm{min}}|\leq |\lambda_{\mathrm{min}}|$ and $|\eta_{\mathrm{max}}|\geq |\lambda_{\mathrm{max}}|$.

It follows that unless your window is a constant magnitude window, the magnitudes of the eigenvalues of the composite operator $FW$ are not constant and cannot be rescaled to 1. Therefore in this case you can always find two different input signals whose magnitudes are changed differently, excluding the possibility of a single normalisation factor.

Or like I commented above: The windowed Fourier transform is not (in general) unitary.

However, you might want to consider if the question was even a sensible one to start with. For all practical purposes you want to preserve the energy of the windowed signal. And this is still guaranteed by the unitarity of the Fourier transform.

$\endgroup$
  • 1
    $\begingroup$ I really like the mathematical yet intuitive approach you made. +1. $\endgroup$ – Tendero Mar 6 '18 at 19:46
  • $\begingroup$ For all practical purposes you want to preserve the energy of the windowed signal. I agree, that's one way of looking at it. Still, many applications make an attempt at computing the energies present in the original (unwindowed) signal. The technique I described is not my own (it's in the article Jojek linked to holometer.fnal.gov/GH_FFT.pdf) $\endgroup$ – dsp_user Mar 6 '18 at 20:17
  • $\begingroup$ @dsp_user I honestly don't believe you understand the question and how Jojek's answer relates to it. Your answer remains wrong for this question, no matter how much you want it to be correct. $\endgroup$ – Jazzmaniac Mar 6 '18 at 20:19
  • $\begingroup$ I provided a link to the pdf article for everyone to see. I think you're misinterpreting what I meant. $\endgroup$ – dsp_user Mar 6 '18 at 20:22
  • $\begingroup$ @dsp_user There is not much to misinterpret. You obviously mean it as an answer to this question. And it isn't. It's the answer to different (sensible) questions you can ask. $\endgroup$ – Jazzmaniac Mar 6 '18 at 20:23
3
$\begingroup$

[EDIT: 20180307, expanded some details] Globaly no, windowing does NOT affect Parseval's theorem (because theorems are only affected, more precisely not applicable, when their hypotheses are not met), in the sense that the equality in energy relates a signal (windowed $x_w$ or not, $x$, resp.) and its Fourier transform (from a windowed signal, $X_w$ or not $X$, resp.). In mathematical terms, you have $\|x\|^2 = \|X\|^2$ and $\|x_w\|^2 = \|X_w\|^2$, but in general, $\|x\|^2 \neq \|X_w\|^2$. The question is: can we correct that, simply? Globaly no again, but approximately so, yes.

First, let us take a little step back. Parseval's theorem, sometimes called Rayleigh's identity or energy theorem, is part of a more generic theorem, the Plancherel theorem, itself related to the Bessel inequality: in an Hilbert space, if $e_{1},e_{2},\ldots$ is an orthonormal sequence, and $\left\langle \cdot,\cdot\right\rangle$ is the scalar product, then:

$$\sum _{k=1}^{\infty }\left\vert \left\langle x,e_{k}\right\rangle \right\vert ^{2}\leq \left\Vert x\right\Vert ^{2}\,.$$

The main message here is that, even for orthogonal transforms, energy preservation is not granted.

There are many flavors of Fourier transforms: continuous, discrete time, discrete, Fourier series, and many others. For the first ones it happens that, their properties entail that, for (many) Fourier transforms, the Bessel equality is indeed an equality (with mild hypotheses on the function).

This is quite important in applications, and their associated algorithms. In simpler words: whenever one removes stuff from the Fourier domain (filtering, deconvolution, restoration, etc.), s(h)e removes the same stuff in the original domain (time, space). Energy conservation, even approximately, plays a great role in algorithmic stability, inverses, etc.

Coming back to windows, their uses in signal processing are legion: apodization data, limiting artifacts, providing approximate stationarity to non-stationary signals, reducing leakage, yielding running/short-time transforms, allowing parallel processing, etc.

So, the present question is, for is the window for, and what are its effects? The span is large. A uniform window is likely to have little effect. A very concentrated window would have a tremendous influence, because it will only concentrate on some signal samples, important of not.

Simply put: take a uniform window, it does not affect samples, and Parseval-Plancherel is preserved. take a null window except on one sample, the energy could vary a lot, depending on where the window is located. At one extremity, if you applied a zero window, the signal would be zero, with zero energy, so you does retain the energy of the original signal. Somehow, the window effect is related to both the window shape, and signal's properties.

Now, when one is sliding windows, like in time-frequency/multirate filter bank processing, each sample at the end sees the window passing, and get an equal share of amplitude. For properly chose windows (with unit energy), correction factors is more easily applied.

For a single signal frame, the outcome will be signal/window dépendant a lot. Let us concentrate on discrete signals $x$ and a given window $w$.

For the FFT side, the window is uniform, and the relationship between the signal energy and the FFT version is linear, as shown in the following graph:

clear;close all;
nSample = 16^2;
nRealization =16^2 ;
nRandWindow = 2^6;

matrixResultSignalEnergy = zeros(nRealization,1);
matrixResultFourierEnergy = zeros(nRealization,1);
matrixResultFourierEnergyRandWindow = zeros(nRealization,1);

for iRealization = 1:nRealization
    data = randn(nSample,1);
    dataFFT = fft(data);
    matrixResultSignalEnergy(iRealization) = norm(data).^2;
    matrixResultFourierEnergy(iRealization) = norm(abs(dataFFT)).^2/nSample;
    for  iRandWindow = 1:nRandWindow
        randWindow = rand(nSample,1);
        randWindow = randWindow/norm(randWindow);
        matrixResultFourierEnergyRandWindow(iRealization,iRandWindow) = norm(abs(fft(data.*randWindow))).^2/nSample;
    end
end

figure(1);clf;hold on
plot(matrixResultSignalEnergy,matrixResultFourierEnergy,'x');
xlabel('Signal energy');
ylabel('Fourier energy');
grid on;axis tight



figure(2);clf;hold on
% errorbar(matrixResultSignalEnergy,(matrixResultFourierEnergy),std(matrixResultFourierEnergy));
errorbar(matrixResultSignalEnergy,mean(matrixResultFourierEnergyRandWindow,2)*nSample,std(matrixResultFourierEnergyRandWindow')/sqrt(nSample),'.');
grid on;axis tight
    grid on;axis tight

Details follow. For a uniform window, the dependence is strictly linear:

enter image description here

So, no big deal. But for non uniform windows, what can happen? For a Gaussian random data, and a window drawn from a positive random, unit energy uniform distribution, you get:

enter image description here

So, the slightly dispersed cloud tells you that there is not single factor you can apply in general to preserve energy with a window, although the relation is close to linear.

In practice, you could employ a more flat window, which reduces the portion of samples (at the edge) that are distorted with the window. Raided cosine windows can be useful in this option.

From an alternative side, DSP people have developed techniques to compensate the amplitude or energy variations, as described in:

$\endgroup$
  • 1
    $\begingroup$ Nice answer, especially because you explained it from both a theoretical and practical aspect. $\endgroup$ – dsp_user Mar 7 '18 at 7:39
  • $\begingroup$ I am glad you appreciate. Should I add details/simulations on intermediate windows between uniform and random? $\endgroup$ – Laurent Duval Mar 7 '18 at 19:41
  • $\begingroup$ Sure, the more information you add, the better (especially for people reading this in the future who might be confused). I'm sorry that we can't all agree what is actually being asked here. I mean, some people (Jazzmaniac primarily) view this whole question from a purely theoretical point of view whereas I'm mostly concerned with how to compute normalized magnitude values (and the OP explicitly asked for that). That's why I think that your answer is important because it tries to reconcile both views. $\endgroup$ – dsp_user Mar 7 '18 at 20:32
1
$\begingroup$

Tendero's answer clearly explains why Parseval's theorem is still valid so I'll just address the question of getting the correct magnitudes (normalization). First off, it should be emphasized that windowing always results in energy loss (except in the case of a rectangular window). In order to compute the correct magnitudes we must compensate for this energy loss. Fortunately, this is rather easy to do.

The first thing to do is to compute the two sums, s1 and s2 respectively, which will be different for every window type. The actual sums will also depend on the FFT/window length (len). These two sums will be used later for scaling the magnitudes (normalization).

//c++ code (using a Hamming window as an exmaple)

float *CreateHammingWindow(unsigned len, float &s1, float &s2){
 s1 = 1;
 s2 = 1;
 float *pWnd = new float[len];
 unsigned lenHalf = len/2;
 //set
 pWnd[0] = 0;
 pWnd[lenHalf] = 1;
 for(int i=1,j=len-1; i<lenHalf; i++,j--){
  pWnd[i] = 0.54 - 0.46 * cos(PI*2 * i/len);//here we just calculate values for a Hamming window
  pWnd[j] = pWnd[i];
  s1 +=  (2 * pWnd[i]);//multiply by 2 because we're using symmetric windows
  s2 +=  (2 * (pWnd[i] * pWnd[i])) ;
 }
 return pWnd; //window data
}

Now, in order to compute a normalized power spectrum (for every bin), we can use the following equation

for(int i=0;iN/2;i++)
    Ps[i] = 2*((FFT_real(i)*FFT_real(i) + FFT_imaginary(i)*FFT_imaginary(i))/(s1*s1));

In order to get magnitudes in decibels, you can use

Psdb = 10*log(Ps[i]) //power spectrum in dB

or

Psdb = 10*log(sqrt(Ps[i])) //linear power spectrum in dB

Note that you may still have to multiply by a constant (C) to get the correct values (could be because of overlapping etc)

Psdb = Psdb * C;

Similarly, s2 is used to compute PSD (power spectral density).

$\endgroup$
  • $\begingroup$ This is very misleading. A windowed Fourier transform is not unitary. An in general there is no way to fix that by scaling. $\endgroup$ – Jazzmaniac Mar 6 '18 at 15:23
  • $\begingroup$ I've found it right here on dsp.stackexchange. Please have a look at Jojek's answer and the article he linked to (dsp.stackexchange.com/questions/32076/…). The code I posted is something I'm using for my own project. So please, make sure you check some facts before you downvote an answer. $\endgroup$ – dsp_user Mar 6 '18 at 17:57
  • $\begingroup$ I stand by what I said. If you think otherwise please show a mathematical proof. Specifically, Jojek's answer only holds for the special case of a stationary signal. $\endgroup$ – Jazzmaniac Mar 6 '18 at 18:18
  • $\begingroup$ Do as you please. Everyone's entitled to their opinion. $\endgroup$ – dsp_user Mar 6 '18 at 18:20
  • 1
    $\begingroup$ the problem @dsp_user, is that you're wrong. $\endgroup$ – robert bristow-johnson Mar 6 '18 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.