FFT code in Fortran

Hi,

I would like to use the Fortran code pasted below in Matlab. Is there an equivalent function in ML for that code? Does fft do the job? Can you help me to "translate" it into ML?

C Subroutine FFT, Cooley-Tukey radix-2, DIF (decimation-in-frequency)
C FFT program written by C. S. Burrus, Rice University, Sept. 1983.
C Complex input data in arrays X (real part) and Y (imaginary part).
      SUBROUTINE FFT(X,Y,N,M)
      REAL X(1),Y(1)
C Main FFT loops
      N2=N
       DO 10 K=1,M
         N1=N2
         N2=N2/2
         E=6.28318531/N1
         A=0.
         DO 20 J=1,N2
            C=COS(A)
            S=SIN(A)
            A=E*J
            DO 30 I=J,N,N1
                L=I+N2
               XT=X(I)-X(L)
               X(I)=X(I)+X(L)
               YT=Y(I)-Y(L)
               Y(I)=Y(I)+Y(L)
               X(L)=C*XT+S*YT
               Y(L)=C*YT-S*XT
   30       CONTINUE
    20    CONTINUE
   10 CONTINUE
C Digit reverse counter
  100 J=1
      N1=N-1
      DO 104 I=1,N1
         IF (I.GE.J) GO TO 101
         XT=X(J)
         X(J)=X(I)
         X(I)=XT
         XT=Y(J)
          Y(J)=Y(I)
         Y(I)=XT
  101    K=N/2
  102    IF (K.GE.J) GO TO 103
         J=J-K
         K=K/2
         GO TO 102
  103    J=J+K
  104 CONTINUE
      RETURN
      END

Thank you for ur answers!

Answer 1

[Image Analyst]

Example of a 2D FFT:

% 2D FFT Demo
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;

% Change the current folder to the folder of this m-file.
if(~isdeployed)
  cd(fileparts(which(mfilename)));
end
  
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

% Read in a standard MATLAB gray scale demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'cameraman.tif';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
% Read in image.
grayImage = imread('cameraman.tif');
[rows columns numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
  grayImage = rgb2gray(grayImage);
end

% Display original grayscale image.
subplot(2, 2, 1);
imshow(grayImage)
title('Original Gray Scale Image', 'FontSize', fontSize)

% Perform 2D FFTs
fftOriginal = fft2(double(grayImage));
shiftedFFT = fftshift(fftOriginal);
subplot(2, 2, 2);
imshow(real(shiftedFFT));
title('Real Part of Spectrum', 'FontSize', fontSize)
subplot(2, 2, 3);
imshow(imag(shiftedFFT));
title('Imaginary Part of Spectrum', 'FontSize', fontSize)

% Display magnitude and phase of 2D FFTs
subplot(2, 2, 4);
imshow(log(abs(shiftedFFT)),[]);
colormap gray
title('Log Magnitude of Spectrum', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);

% Now convolve with a 2D rect function.
figure;
rectWidth = 10;
rectHeight = 5;
kernel = ones(rectHeight, rectWidth) / (rectHeight * rectWidth);
% Display it
subplot(2, 2, 1);
k = padarray(kernel, [3, 3]); % Just for display.
imshow(k, []);
axis on;
title('Kernel', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);

% Convolve kernel (box filter) with the image
filteredImage = conv2(double(grayImage), kernel, 'same');
% Display filtered image.
subplot(2, 2, 2);
imshow(filteredImage,[]);
title('Filtered Image', 'FontSize', fontSize)

% Perform 2D FFT on the filtered image to see its spectrum.
% We expect to see a sinc multiplication effect.
% It should look like the original but with a sinc pattern overlaid on it.
fftFiltered = fft2(double(filteredImage));
shiftedFFT = fftshift(fftFiltered);

% Display magnitude of the 2D FFT of the filtered image.
subplot(2, 2, 3);
imshow(log(abs(shiftedFFT)),[]);
colormap gray
title('Log Magnitude of Spectrum - Note sinc multiplication', 'FontSize', fontSize)

Answer 2

[Rick Rosson]

Please try the following:

   Fs = 48000;
   dt = 1/Fs;
   t = (0:dt:0.25-dt)';

   N = size(t,1);
   dF = Fs/N;
   f = (-Fs/2:dF:Fs/2-dF)' + mod(N,2)*dF/2;

   Fc = 500;
   x = cos(2*pi*Fc*t);
   y = sin(2*pi*Fc*t);
   z = x + 1j*y;

   Z = fftshift(fft(z))/N;

   figure;

   subplot(2,1,1);
   plot(t,x,t,y);

   subplot(2,1,2);
   plot(f,abs(Z));

Answer 3

[Wayne King]

Why not just use fft() in MATLAB?

It is a built-in function so I'm sure it would be just as efficient (probably more) than creating a MEX file from the Fortran source. If you really want to learn to create MEX files from Fortran source see

http://www.mathworks.com/help/matlab/creating-c-c-and-fortran-programs-to-be-callable-from-matlab-mex-files.html

[Rok mr]

Thank you for ur answer!

Can u help me with an example how to use fft function instead the Fortran subroutine. X are real parts of complex vector and Y are imaginary parts. N ... number of samples ; N = 2^M

Does this give me the same result as the Fortran subroutine: Xr = fft(X+1i*Y);

I have some problems with the result scaling in my program and I would like to check if it comes from the fft.

how is it if I want to use the Fortran subroutine for inverse fft? Is the difference only that I have to multiply the result by 2/N?

[Image Analyst]

Apparently you didn't look in the help. There is an example right in there. Plus, Rick gave you another example below. Of course more examples are always better so maybe Wayne will give one. If you want a 2D fft example, see my answer.

[Rick Rosson]

   >> doc fft
   >> doc ifft
   >> doc fftshift
   >> doc abs
   >> doc log10
   >> doc complex

Answer 4

[Rick Rosson]

First, you do not need to use loops in MATLAB to compute the values of vectors and matrices. Second, if you want to compute the inverse Fourier transform, you must provide the two-sided spectrum, from -OMAX to OMAX (or from 0 to 2*OMAX), not the one-sided spectrum from 0 to OMAX.

Please try the following:

   % Time increment (seconds per sample):
   dt = pi/OMAX;

   % Time domain (seconds):
   t = dt*(0:N-1)';

   % Sampling rate (samples per second):
   Fs = 1/dt;

   % Frequency increment (hertz per sample):
   dF = Fs/N;

   % Frequency domain (hertz):
   f = (-Fs/2:dF:Fs/2-dF)' + mod(N,2)*dF/2;

   % Complex angular frequency (radians per second):
   omega = 2*pi*f; 
   s = 1j*omega;

   % Complex impedance, frequency response ???
   Z = C*s .* (L*s + R) ;
   ...
   ...