fft and inverse fft error size and applying impulse response matlab

Question

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.

Answer 1

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

Answer 2

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

Answer 3

[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.