1
$\begingroup$

I have a signal called noise. I am using the following command to calculate the mean squared error after fft and inverse fft process:

 sum((noise-real(ifft(fft(noise,2048),num_of_samples))).^2)/num_of_samples

Where num_of_samples=2819519 is the number of samples in the original signal. I am getting an error of $\sim10^{-8}$. Is that reasonable and due to approximation or does that indicate something is not converted correctly?

My final goal is to multiply the signal in the frequency domain with a $1052\times 4$ given impulse response. This represents a room with 4 microphones.

$\endgroup$
  • $\begingroup$ is the num_of_samples = 2048? $\endgroup$ – Fat32 Jun 16 at 13:38
  • $\begingroup$ No. It is the same number of samples as the original signal. $\endgroup$ – havakok Jun 16 at 13:55
  • $\begingroup$ what is that number ? $\endgroup$ – Fat32 Jun 16 at 13:57
  • $\begingroup$ 2819519. Does it have any impact? $\endgroup$ – havakok Jun 16 at 13:58
  • $\begingroup$ yes... you are taking 2048 point FFT of 2819519 point sequence and trying to convert it back ? $\endgroup$ – Fat32 Jun 16 at 13:59
1
$\begingroup$

I don't know your overall goal but if you want to test the MSE associated with FFT / IFFT then you should perform the following

num_of_samples = 2819519 ;
fft_len         = 2^nextpow2(2819519) ;

X = fft(signal, fft_len);
y = real( ifft(X, fft_len) );
mse = sum( (signal - y(1:num_of_samples) ).^2 )/num_of_samples
$\endgroup$
  • $\begingroup$ havakok, so what's your MSE now ? $\endgroup$ – Fat32 Jun 16 at 17:09
  • $\begingroup$ hay Fat, is that you in the pic? $\endgroup$ – robert bristow-johnson Jun 17 at 6:24
  • $\begingroup$ MSE is good now, though I am still trying to figure out how to multiply this in the frequency domain with the given size. $\endgroup$ – havakok Jun 17 at 6:32
  • 1
    $\begingroup$ @robertbristow-johnson probably so ;-) $\endgroup$ – Fat32 Jun 17 at 8:35
1
$\begingroup$

Try the snippet below.

In double precision the relative error should be around -300 dB. The absolute error depends on the scale of your signal.

From your post it sounds you want to filter a long signal with a bunch of much shorter impulse responses. The best way to do this would be "overlap add", which is implemented by MATLAB's fftfilt function.

nx = 2.^(4:16);        % gets some FFT sizes
for i=1:length(nx)   
  x  = randn(nx(i),1); % create noise
  y= ifft(fft(x));     % FFT and IFFT
  d = x-y;             % difference 
  % calculate relative difference in dB
  err = 10*log10(sum(d.^2)./sum(x.^2));
  fprintf('FFT Length = %5d, Error = %7.2fdB\n',nx(i),err);
end
$\endgroup$
0
$\begingroup$

[General considerations before trying to answer] Several artifacts may happen when one implements a "theoretically" inversible transform with finite-size computing (integer, float, double float, etc.). In Floating point error mitigation, you can read:

Floating-point error arises because real numbers cannot, in general, be accurately represented in a fixed space. By definition, floating-point error cannot be eliminated, and, at best, can only be managed.

H. M. Sierra noted in his 1956 patent "Floating Decimal Point Arithmetic Control Means for Calculator":

"Thus under some conditions, the major portion of the significant data digits may lie beyond the capacity of the registers. Therefore, the result obtained may have little meaning if not totally erroneous."

The classical caveat is often stated as $1/3=0.3333333\ldots$ but $3\times 0.3333333\ldots \neq 1$. The classical group or ring rules of arithmetics like:

  • inverse: $a.a^{-1}=1$ (see above)
  • associativity/commutativity: $(a+a)+(b+b) = (a+b)+(a+b)$

are not always satisfied on computers. For the second one, this may happen in multi-core or parallel computing, depending on which operation is performed on which processor. This becomes even more imprecise as the number of operands (the values you operate on) increases. This could even become non-deterministic, as discussed in Can floating point error (in FFTW3) cause non-deterministic behavior?. Even with uint8 data, you can get unexpected results, due to rounding (What do the 1D filters represent when using imfilter?).

[Back to the question] Padding the FFT or truncating its inverse are not linear opearations. However, I will forget them for now. But sum, FFT, IFFT, real-part or division by the number of samples are all linear. If I am not wrong in taking a step back, what you are trying to estimate for a real signal $s(t)$, $\mathcal{F}$ denoting Fourier, is:

$$ E_{T_1}^{T_2}(s) = \frac{1}{T_2-T_1}\int_{T_1}^{T_2}\left(s(t) - \Re\left(\mathcal{F}^{-1} \left(\mathcal{F}(s(t)\right)\right)\right)^2 dt$$

which should be zero for a real signal. So, first hypothesis: your signal is not real. And the first action should be to test this.

If this is settled, you can easily infer that

  • if you multiply $s$ by $\lambda$, you get $E_{T_1}^{T_2}(\lambda s) = |\lambda|^2 E_{T_1}^{T_2}( s) $
  • if the diffference $s(t) - \Re\left(\mathcal{F}^{-1} \left(\mathcal{F}(s(t)\right)\right)$ is "somehow second order stationnary", the energy is somehow linear with the time-span.

If you convert it to discrete samples, you get your Matlab formula, which is (roughly):

  • quadratic wth the scale (amplitude) or our signal,
  • inversely linear with the integrated length.

So, in a word, you are considering a somehow absolute error of $∼10^{−8}$, which is thus VERY difficult to interpret, as it related to $S^2/N$, if $S$ is the discrete signal's amplitude, and $N$ its length, as discussed by @Hilmar. Increase the length, the value will decrease.

If you expect the residue to be stationary, a relative error is much more informative. On my side, I usually relate to the number of bits data is stored on and computations are performed.

$\endgroup$

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.