FFT code in Fortran
31 次查看(过去 30 天)
显示 更早的评论
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!
0 个评论
回答(4 个)
Image Analyst
2012-12-25
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)
0 个评论
Rick Rosson
2012-12-26
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));
0 个评论
Wayne King
2012-12-23
编辑:Wayne King
2012-12-23
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
3 个评论
Image Analyst
2012-12-25
编辑:Image Analyst
2012-12-25
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
2012-12-26
>> doc fft
>> doc ifft
>> doc fftshift
>> doc abs
>> doc log10
>> doc complex
Rick Rosson
2012-12-24
编辑:Rick Rosson
2012-12-24
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) ;
...
...
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!