Re: FFT spectrum analyzer question

[Bolin Hsu]

Macy,

Thanks for your reply. Strangely I get 0.15578 instead of 1/2 from the following:

Ss=192000;
dt=1/Ss;
N=9600;
t=[0:N-1]*dt;
sig=sin(2*pi()*t);
sig_f=fft(sig)/N;
max(abs(sig_f))

Anyway, I was studying Bartlett's method for power spectrum estimation. That's why I fixed a run length, split the signal into multiple sections, and apply fft on these sections. If I apply the same idea on your sample code, I get the following:

Ss=192000;
dt=1/Ss;

for N = [192000 96000 48000 24000 12000 6000]
  max_sum = 0;
  sections = 0;
  for offset = [1:N:Ss]
    t=[offset:offset+N-1]*dt;
    sig=sin(2*pi*t);
    sig_f=fft(sig)/N;
    max_sum = max_sum + max(abs(sig_f));
    sections = sections + 1;
    max_avg = max_sum / sections;
    printf('[%6d %6d] max: %f avg: %f\n', N, i, max(abs(sig_f)), max_avg)
  endfor
  printf('\n')
endfor

The test shows that max(abs(sig_f)) is 1/2 only when N is equal to Ss. And the average of the max over all sections doesn't converge to 1/2, either. Is this inherent in Bartlett's method or caused by bugs in my implementation?

Thanks,

Bolin

On Wed, May 14, 2014 at 8:31 AM, Macy <address@hidden> wrote:

lot's of 'errors'

given signal length as len [I prefer labeling run length as N]
sig_f=fft(sig)/len;
you do that you'll get closer to constant.

next be sure to make the length 'correct' AND match sample lengths, must be ineteger related else you will have to use 'windowing' as partial values get shoved into adjacent bins

your frequency, sampling rate, and packet lengths 'fight' each other.

also, for time use one less than the length you set
to understand try this:
Ss=192000;
dt=1/Ss;
N=9600;
t=[0:N-1]*dt;
sig=sin(2*pi()*t);
sig_f=fft(sig)/N;
max(abs(sig_f))
will always come up 1/2
since the energy gets split during your fft.
if you plot
plot(20*log10(abs(sig_f)));
you'll notice two halves to your spectrum which can get a bit confusing.
just remember that your time waveform is ALWAYS real so the two havlesof the resulting spectrum will ALWAYS be related as conjugates, BUT the right half is inverted, coming down. see what i mean about confusing?

now do again using
t=[0:N]; and the plots will suck.


I wrote a way to bounce back and forth between time and frequency domains:
two custom functions: time2freq and freq2time
[sig_f,ford]=time2freq(sig);
so when you plot
plot(ford/dt,20*log10(abs(sig_f)));
you get an ACTUAL spectrum plot.

if you want log-log
plot(log10(ford(2:end)/dt/fc),20*log10(abs(sig_f(2:end))));
why not index 1? because that is DC and it plays havoc with your log functions.

The advanatage of this form is that you can do all kinds of spectral manipulations like run your signal through a filter by multiplying the two spectrums, then convert back to time and see the effect on the time waveform.

sig_t=freq2time(sig_f.*filterspectrum);
and get waveform back.

ps: always use even length runs, gets really confusing if you don't.






--- address@hidden wrote:

From: Bolin Hsu <address@hidden>
To: address@hidden
Subject: FFT spectrum analyzer question
Date: Tue, 13 May 2014 17:09:59 -0700

I saw the thread about FFT spectrum analyzer at
http://octave.1599824.n4.nabble.com/FFT-Spectrum-Analyzer-td4630647.html,
and gave it a try. From http://en.wikipedia.org/wiki/Bartlett%27s_method, I
got the idea that the energy is abs(fft(x))^2/M. So I ended up with the
program below. The signal was sine wave with peak amplitude 1, so I thought
I would get the RMS amplitude 1/sqrt(2). But the result I got was much
smaller. And the result varies with different FFT size. What am I missing
in this program?

Fs=44100;
w=14000;
t=0:1/Fs:1;
x=sin(2*pi*w*t);
len=length(x);
for N = [128 256 512 1024 2048 4096 len]
  for i = [1:N:len]
    if (i+N-1 <= len)
      xdft=fft(x(i:i+N-1));
      xdft=xdft(1:N/2+1);
      psdx=(1/(Fs*N)).*abs(xdft).^2;
      psdx(2:end-1)=2*psdx(2:end-1);
      freq=0:Fs/N:Fs/2;
      printf('[%5d, %5d] %f\n', N, i, max(psdx))
    endif
  endfor
endfor


Thanks,
Bolin

_______________________________________________
Help-octave mailing list
address@hidden
https://lists.gnu.org/mailman/listinfo/help-octave