Discussion Forums

Home › Forums › Signal Processing Algorithms

Octave FFT Average function

7 replies [Last post]
Thu, 2010-07-22 17:42
af6kd
Offline
Joined: 2010-04-25
Posts: 23

Hello,

Below is a script that I have been using to create an averaged and windowed FFT:

function ff = fft_avg(x, M)

  if (nargin == 1)
     M = length(x);
  endif;

N = length(x);

win = blackman(M);
%win = hanning(M);
%win = hamming(M);

%% calculate average FFT, with window size M

ff=zeros(M,1);
for ci=1:N/M,
 ff=ff+abs(fft(x((M*(ci-1)+1:M*ci)).*win));
end
ff=ff/N*M;

This has been cobbled together from multiple scripts, and if anyone has any input, it would be appreciated.

Dave

Top
Thu, 2010-07-22 23:13
#1
gerryharp
Offline
Joined: 2010-05-15
Posts: 364
window (apodization) is good

Hi Dave
Thanks for contributing this code. We were planning to put in windowing in front of our power spectra, next.
 
Another way to improve the powerspectrum example is to use, say, a Hanning window and then step through the data only 1/2 of an FFT_length at a time. As you know, window functions reduce boundary artifacts in your FFT, but they effectively throw away 1/2 of the data (SNR drops by root(2), your mileage varies depending on which window you use). The correct way to handle this deficiency is to perform 2x as many FFT's, stepping only one half as far each time. By reducing the step size, you make sure no data is lost and recover the correct total power at the end.
 
Want to modify the Ocave code and post it?
 
Gerry

Top
Thu, 2010-07-22 23:28
#2
af6kd
Offline
Joined: 2010-04-25
Posts: 23
Good idea, I'll try to get to

Good idea, I'll try to get to that soon....    I'll see if I can't work it into your Power Sepctrum code too.

While I like the new forum format, is this the best way to post code?  Sometimes it is hard to find that snippet of code that lurks deep within the forum...

Dave

Top
Fri, 2010-07-23 05:45
#3
gerryharp
Offline
Joined: 2010-05-15
Posts: 364
posting code

Hi Dave
 
Clearly this is not the best way to post code. I would like to have the ability to attach documents to posts. Let me ask you and the community: Would document attachements be the best way to post code snippets (or other things, like graphics)?
 
Gerry

Top
Fri, 2010-07-23 05:56
#4
Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Code that is not suitable for

Code that is not suitable for being posted in a source repository I prefer to get in some kind of plain text form (I've used pastebin.com here on the forum a couple of times). Graphics I prefer to get as PNG or JPG, so they can be accessed directly in the browser (imageshack.us and similar services will host them).

Something like this would be ideal, but it may not be easy to set up.

Top
Fri, 2010-07-23 16:49
#5
af6kd
Offline
Joined: 2010-04-25
Posts: 23
If I look at the other forums

If I look at the other forums I read, attachments are often used.   It would be nice if graphics/plots would appear without having to click a link.

I guess for small code snippets, the forum is OK (that is common on other Matlab/octave forums).  

How about a GIT repository?  

Top
Fri, 2010-07-23 17:03
#6
Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
 I agree. Not only can code

 I agree. Not only can code in the repository easily be linked to (with syntax highlighting and everything), it can also easily be downloaded to one's own workstation, even when the code spans several source files.

Top
Sat, 2010-07-24 01:38
#7
af6kd
Offline
Joined: 2010-04-25
Posts: 23
Updated FFT averaging function

Here is an updated FFT averaging function, rewritten and made easier to understand. 

Again, feedback is appreciated!  I *think* it is working correctly now, with the 1/2 scale windowing.  :)

Dave

function y = fft_avg2(x, win_length)

  if (nargin == 1)
    win_length = 4096;
  endif;

step = win_length / 2;

fft_win = hanning(win_length)';

y = zeros(1,win_length);

start = 1;
stop = win_length;

while (stop < length(x) - step)
 y = y + abs(fft(x(start:stop).*fft_win));
 start = start + step;
 stop = stop + step;
end
 
y=y/(length(x)/step);

Top