help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

strange error inside loop in octave 3.6.1 (mingw)


From: david moloney
Subject: strange error inside loop in octave 3.6.1 (mingw)
Date: Wed, 15 Aug 2012 10:43:20 +0100

The radix4 fft code below works fine is i set power of 4 (xp) values case-by-case.

$ octave fft4.m
GNU Octave, version 3.6.1
Copyright (C) 2012 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.

Octave was configured for "i686-pc-mingw32".

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html

Read http://www.octave.org/bugs.html to learn how to submit bug reports.

For information about changes from previous versions, type `news'.

 - Use `pkg list' to see a list of installed packages.
 - MSYS shell available (c:\Program Files (x86)\Octave-3.6.1\msys).
 - Graphics backend: qt.

ans = 1.4198e-015

However if I uncomment the loop code I get the following error

$ octave fft4.m
GNU Octave, version 3.6.1
Copyright (C) 2012 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.

Octave was configured for "i686-pc-mingw32".

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html

Read http://www.octave.org/bugs.html to learn how to submit bug reports.

For information about changes from previous versions, type `news'.

 - Use `pkg list' to see a list of installed packages.
 - MSYS shell available (c:\Program Files (x86)\Octave-3.6.1\msys).
 - Graphics backend: qt.

error: `stage' undefined near line 48 column 68
error: evaluating argument list element number 1
error: evaluating argument list element number 2
error: called from:
error:   r4fftN at line 48, column 22
error:   c:\Users\david\Documents\Visual Studio 2010\Projects\mv_fft\fft4.m at l
ine 80, column 7

the "error" refers to a line the in fft function code which otherwise works correctly when xp is not set by a loop ... very strange.


-------------------------------------------------- fft4.m octave/matlab code ---------------------------------------------------------


function Z = radix4bfly(x,segment,stageFlag,W)
  % For the last stage of a radix-4 FFT all the ABCD multiplers are 1.
  % Use the stageFlag variable to indicate the last stage
  % stageFlag = 0 indicates last FFT stage, set to 1 otherwise

  % Initialize variables and scale to 1/4
  a=x(1)*.25;
  b=x(2)*.25;
  c=x(3)*.25;
  d=x(4)*.25;

  % Radix-4 Algorithm
  A=a+b+c+d;
  B=(a-b+c-d)*W(2*segment*stageFlag + 1);
  C=(a-b*j-c+d*j)*W(segment*stageFlag + 1);
  D=(a+b*j-c-d*j)*W(3*segment*stageFlag + 1);

  % assemble output
  Z = [A B C D];
end % radix4bfly()


% radix-4 DIF FFT, input signal must be floating point, real or complex
%
function S = r4fftN(s)

  % Initialize variables and signals: length of input signal is a power of 4
  N = length(s);
  M = log2(N)/2;

  % Initialize variables for floating point sim
  W=exp(-j*2*pi*(0:N-1)/N);
  S = complex(zeros(1,N));
  sTemp = complex(zeros(1,N));

  % FFT algorithm
  % Calculate butterflies for first M-1 stages
  sTemp = s;
  for stage = 0:M-2
    for n=1:N/4
      S((1:4)+(n-1)*4) = radix4bfly(sTemp(n:N/4:end), floor((n-1)/(4^stage)) *(4
^stage), 1, W);
    end
    sTemp = S;
  end

  % Calculate butterflies for last stage
  for n=1:N/4
    S((1:4)+(n-1)*4) = radix4bfly(sTemp(n:N/4:end), floor((n-1)/(4^stage)) * (4^
stage), 0, W);
  end
  sTemp = S;

  % Rescale the final output
  S = S*N;
end % r4fftN(s)


% test FFT code
%
xp = 2;
%for xp=1:8
  N = 4^xp; % must be power of: 4 16 64 256 1024 4086 ....
  x = 2*pi/N * (0:N-1);
  x = cos(x);
  Y_ref = fft(x);
  Y = r4fftN(x);
  Y = digitrevorder(Y,4);
  %Y = bitrevorder(Y,4);
  abs(sum(Y_ref-Y)) % compare fft4 to built-in fft
%end

reply via email to

[Prev in Thread] Current Thread [Next in Thread]