Help with simple FFT problem

Hi all, Im taking my first steps in FTT processing, and I have a simple doubt that I cant figure out. Im trying a 1.41V sine signal (ie. 1V rms)@1khz and then doing the FFT, i expected to see a 94dBpeak at 1khz but the actual value is incorrect when i plot the figure. Can anybody help me please?? thanks!!

Here's the code:

t=0:1/44100:2;
fs=44100;

VarName1=1.41.*sin(2.*pi().*1000.*t);
%1.41pa=1pa rms=94dB

Pref=20e-6;

N=length(VarName1);
RMS2=sqrt(mean(VarName1.^2))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DbTOTAL=20*log10(RMS2/Pref)

y=fft(VarName1); 

%why the peak at 1000hz is not 94dB?
db = 20*log10((2/N).*(abs(y)/Pref).*(RMS2));

f = (0:length(t)-1)*fs/length(t);
figure
plot(f,db);

Answer 1

[Matt J]

The fft needs to be normalized according to the sampling to make y of the same scale as the continous Fourier transform

y=fft(VarName1)/fs;

There may still be some inaccuracy due to aliasing and truncation of the sinusoid.

[Matt J]

Also, your frequency samples don't precisely coincide with the integers like 1000 Hz. The frequency sampling interval is very close to 0.5, but not exactly equal to it,

    >> fs/N-.5

ans =

    -5.6689e-06

[Matt J]

You need to be more careful about setting up your sampling. You should set up both your frequency and your time sampling axes FIRST before doing any other DSP

            N=88200;
            df=0.5;    %frequency sampling interval
            dt=1/N/df; %time sampling interval         
            fs=1/dt;

            t=(0:N-1)*dt;
            f=(0:N-1)*df;

      %%spectral analysis last

       VarName1=1.41.*sin(2.*pi.*1000.*t);
       y=abs(fft(VarName1))*dt; 

       plot(f,y);

Answer 2

[Youssef Khmou]

You have to change the sampling rate, try the following :

fs=4000;
t=0:1/fs:2-1/fs;

 VarName1=1.41.*sin(2.*pi.*1000.*t);
%1.41pa=1pa rms=94dB
Pref=20e-6;
N=length(VarName1);
RMS2=sqrt(mean(VarName1.^2))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DbTOTAL=20*log10(RMS2/Pref)
N=900; % Choose value for FFT number of points 
y=fft(VarName1,N);
f=((0:N-1)/(N-1))*fs;
figure; plot(f,20*log10(abs(y)));

Answer 3

[Geoff Hayes]

Allyson - there is a thread at Matlab wavread and FFT that seems to have a similar inquiry to yours, or rather the final response from Steve Amphlett includes a solution to the question ...the 1khz tone is a 94 db calibration tone. So then how do i use this information to obtain the right db values?. If you follow through his response, he references your Pref and RMS, but uses them in a slightly different manner.

I'll try and summarize what he is doing, and will also use some of the code from section 2.1 Calbiration of Implementing Loudness Models in MATLAB as it is similar to what he (and you) are doing, especially with respect to variable names (though the latter seems to calculate the RMS on the FFT data)

 % set the sample rate
 fs = 44100;

 % set the time period (note that we subtract 1/fs so that the
 % size of t is 2*44100)
 t=0:1/fs:2-1/fs;

 % get the time-domain data
  x = 1.41.*sin(2.*pi().*1000.*t);

 % set the sound pressure level measure
 SPLmeas = 94;

 % set the reference pressure
 Pref    = 20e-6;

 % calculate the RMS
 RMS     = sqrt(mean(x.^2));

 % calculate the calibration factor
 cal = Pref*(10.^(SPLmeas/20))/RMS;

He then applies this calibration factor to the time-domain data, but we can do so the frequency domain data given that it is just a constant.

 % do the FFT
 N = length(x);
 y = fft(x,N); 

 % apply the calibration factor and scale the FFT to be RMS 
 % (multiply by 2/N to get amplitude and then divide by sqrt(2) 
 % to get RMS)
 yRms = cal*abs(y)*(2/N)/sqrt(2);

 % divide by  20e-6 and 20log10 it to get dB
 dbData = 20*log10(yRms/Pref);

 % plot
 freq = (0:length(t)-1)*fs/length(t);
 figure
 plot(freq,dbData);

The figure shows the 94 dB for the 1000 Hz signal. There are some slight differences between the above and the noted section of the linked pdf, so perhaps others will be able to explain them.