fft(single) is giving inconsistent results

I have found that if I have an array of singles and I want to run it through fft so the transform is performed independently for each column, I get a different answer if I add on additional columns. I haven't been able to reproduce using random inputs, but i did save my inputs out from my workspace so the issue can be reproduced. the input time histories are attached.

I basically have a 10000x20 array of single precision, real numbers. If I add on an additional 10 columns of single precision, real numbers and run it through fft, I get a different result than if I just run the 20 columns through fft. There seems to be some interaction between the columns during the calculation, and the transform should definitely be performed independently for each channel.

I am running on 2012b 64-bit, I see the same result on 2010b 64-bit and 2012b 32-bit but NOT on 2010b 32-bit

load('fft vars.mat');
xxx=[bData mData];
isequal(bData(:,2),xxx(:,2))
f1=fft(xxx);
f2=fft(bData);
f1(591,2)
f2(591,2)

%what?!!

ans =

     1

ans =

 -40.3370 -39.9248i

ans =

 -38.6038 -39.4498i

[Jeremy]

All three of the input variables are single precision, the outputs are all complex single. I do realize that at the frequency I am checking the magnitude is 150dB less than the maximum, so I would expect some loss of resolution, but I would think it should be consistent.

[Matt J]

Attach 'fft vars.mat', so that we can try it, too.

Answer 1

[Matt J]

Well, it is interesting that there are differences when the inputs are supposedly identical. However, the differences are tiny and seem like they must be attributable to floating point round-off issues. After all, the values of abs(f1) and abs(f2) range up to 10^8, so differences between

 -40.3370 -39.9248i

and

 -38.6038 -39.4498i

are well on the fringe of single floating point precision.

[Matt J]

I suspect that the differences might be due to multi-threading. The larger array xxx is partitioned differently among cores when it is processed and so the additions/subtractions involved in the ffts are done in a different order. This would lead to different patterns of floating point errors.

[Jeremy]

yes, I guess I cannot read into it to much. I did direct DFT calculation and it was right between those two answers.