MATLAB: Fft2 function in matlab

MATLAB

Does anyone have idea or good tutorials for manually programming a fft2 function (discrete fast fourier transform)?

The function already exists in MATLAB, I just want to be able to understand how it works

Best Answer

Here is a quick example. I assume a random 2D image where the horizontal axis (columns) represents data collected in space domain and the vertical axis (rows) represents data collected in time domain. The two domains could easily be collected in the same domain. I also chose the number of samples collected in either domain to be different (i.e., a rectangular image), but the number of samples could easily be the same for both dimensions.

Nx = 64; % Number of samples collected along first dimension
Nt = 32; % Number of samples collected along second dimension
dx = 1;  % Distance increment (i.e., Spacing between each column)
dt = .1; % Time increment (i.e., Spacing between each row)
x = 0 : dx : (Nx-1)*dx; % distance
t = 0 : dt : (Nt-1)*dt; % time
data_spacetime = randn(Nt,Nx); % random 2D matrix
Nyq_k = 1/(2*dx); % Nyquist of data in first dimension
Nyq_f = 1/(2*dt); % Nyquist of data in second dimension
dk = 1/(Nx*dx);   % Wavenumber increment
df = 1/(Nt*dt);   % Frequency increment
k = -Nyq_k : dk : Nyq_k-dk; % wavenumber
f = -Nyq_f : df : Nyq_f-df; % frequency
data_wavenumberfrequency = zeros(size(data_spacetime)); % initialize data
% Compute 2D Discrete Fourier Transform
for i1 = 1 : Nx
for j1 = 1 : Nt
for i2 = 1 : Nx
 for j2 = 1 : Nt
  data_wavenumberfrequency(j1,i1) = data_wavenumberfrequency(j1,i1) + ...
   data_spacetime(j2,i2)*exp(-1i*(2*pi)*(k(i1)*x(i2)+f(j1)*t(j2)))*dx*dt;
 end
end
end
end
figure;
subplot(3,1,1);
imagesc(k,f,abs(data_wavenumberfrequency));
colorbar; v = caxis;
title('2D DFT');
fft2result = fftshift(fft2(data_spacetime))*dx*dt;
subplot(3,1,2);
imagesc(k,f,abs(fft2result));
colorbar; caxis(v);
title('FFT2');
subplot(3,1,3);
imagesc(k,f,abs(fft2result-data_wavenumberfrequency));
colorbar; caxis(v);
title('Difference');

Related Solutions

MATLAB: Derivative using FFT properties

Based on the website below:

http://www.cl.cam.ac.uk/teaching/2000/ContMaths/JGD-notes/node11.html

I believe you would do something similar to the example below:

Nx = 64; % Number of samples collected along first dimension

Ny = 32; % Number of samples collected along second dimension

dx = 1;  % x increment (i.e., Spacing between each column)

dy = .1; % y increment (i.e., Spacing between each row)

x = 0 : dx : (Nx-1)*dx; % x-distance

y = 0 : dy : (Ny-1)*dy; % y-distance

data_spacedomain = randn(Ny,Nx); % random 2D matrix

Nyq_kx = 1/(2*dx); % Nyquist of data in first dimension

Nyq_ky = 1/(2*dy); % Nyquist of data in second dimension

dkx = 1/(Nx*dx);   % x-Wavenumber increment

dky = 1/(Ny*dy);   % y-Wavenumber increment

kx = -Nyq_kx : dkx : Nyq_kx-dkx; % x-wavenumber

ky = -Nyq_ky : dky : Nyq_ky-dky; % y-wavenumber

data_wavenumberdomain = zeros(size(data_spacedomain)); % initialize data
% Compute 2D Discrete Fourier Transform
for i1 = 1 : Nx
for j1 = 1 : Ny
for i2 = 1 : Nx
for j2 = 1 : Ny
data_wavenumberdomain(j1,i1) = data_wavenumberdomain(j1,i1) + ...
(2i*pi*ky(j1))*(2i*pi*kx(i1))* ...
(data_spacedomain(j2,i2)*exp(-1i*(2*pi)*(kx(i1)*x(i2)+ky(j1)*y(j2))));
end
end
end
end
differentiated_data = ifft2(ifftshift(data_wavenumberdomain));

I have a program for performing this on a 1D dataset that you might leverage (in your own way) to produce the same effect on 2D datasets - see my answer in the post below:

http://www.mathworks.com/matlabcentral/answers/21700-finding-the-velocity-from-displacement

[EDIT 6/22]

Practically speaking, if you were going to compute the derivative of a 2D image, where your dimensions are either space-space or time-time, then you would need to use fft2, ifftshift, and meshgrid:

Nx = 64; % Number of samples collected along first dimension
Ny = 32; % Number of samples collected along second dimension
dx = 1;  % x increment (i.e., Spacing between each column)
dy = .1; % y increment (i.e., Spacing between each row)
x = 0 : dx : (Nx-1)*dx; % x-distance
y = 0 : dy : (Ny-1)*dy; % y-distance
data_spacedomain = randn(Ny,Nx); % random 2D matrix
Nyq_kx = 1/(2*dx); % Nyquist of data in first dimension
Nyq_ky = 1/(2*dy); % Nyquist of data in second dimension
dkx = 1/(Nx*dx);   % x-Wavenumber increment
dky = 1/(Ny*dy);   % y-Wavenumber increment
kx = -Nyq_kx : dkx : Nyq_kx-dkx; % x-wavenumber
ky = -Nyq_ky : dky : Nyq_ky-dky; % y-wavenumber
% Compute 2D FFT
data_wavenumberdomain = fft2(data_spacedomain);
% Compute grid of wavenumbers
[KX, KY] = meshgrid(ifftshift(kx),ifftshift(ky));
% Compute 2D derivative
data_wavenumberdomain_differentiated = (2i*pi)^2*KX.*KY.*data_wavenumberdomain; 
% Convert back to space domain
data_spacedomain_differentiated = ifft2(data_wavenumberdomain_differentiated );

If you only want to take the partial derivative of the the "image" with respect to one of your dimensions, then all you would do is:

d_MAT_x = 2i*pi*KX.*data_wavenumberdomain;
d_MAT_y = 2i*pi*KY.*data_wavenumberdomain;

Then just transform back into time domain using ifft2

[EDIT 6/27]

1. "Nx" and "Ny" are the number of columns and number of rows, respectively. Thus, you should have:

[KX,FQ]=meshgrid(ifftshift(kx),ifftshift(fq));

2. Your "fq" and "kx" may be defined wrong (I said in an earlier comment that if the number of rows or columns is an odd number then we need to do something different). Also, it seems strange that you would define "fq" in terms of "dx" and "dt"... and "kx" in terms of "dy" and "dx". Lets assume "dt" is the time increment and "dx" is the spatial increment... "Nt" is the number of time samples and "Nx" is the number of spatial samples, then:

Nyq_fq = 1/(2*dt);
Nyq_kx = 1/(2*dx);
dfq = 1/(Nt*dt);
dkx = 1/(Nx*dx);
if mod(Nt,2) == 0 % Nt is even
     fq = -Nyq_fq : dfq : Nyq_fq-dfq;
else            % Nt is odd
     fq = [sort(-1*(dfq:dfq:Nyq_fq)) (0:dfq:Nyq_fq)];
end
if mod(Nx,2) == 0 % Nx is even
     kx = -Nyq_kx : dkx : Nyq_kx-dkx;
else            % Nx is odd
     kx = [sort(-1*(dkx:dkx:Nyq_kx)) (0:dkx:Nyq_kx)];
end

Does this help solve things?

MATLAB: Manual spectrogram creation without using spectrogram command

hello

see below (myspecgram)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
[data,Fs] = audioread('myWAVaudiofile.wav');
channel = 1;
signal = data(:,channel);
samples = length(signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096;    % 
OVERLAP = 0.75;
% spectrogram dB scale
spectrogram_dB_scale = 100;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 2;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
[sg,fsg,tsg] = myspecgram(signal,NFFT,Fs,OVERLAP);  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([' Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function  [sg,freq_vector,time] = myspecgram(signal,nfft,Fs,Overlap)
% [sg,freq_vector,time] = myspecgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
dt = 1/Fs;
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    sg = [];
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        sg = [sg (abs(fft(sw))*4/nfft)];     % X=fft(x.*hanning(N))*4/N; % hanning only 
        time(i) = (start+nfft/2)*dt; % time defined as in the middle of the data buffer
    end
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
    
sg = sg(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Best Answer

Related Solutions

MATLAB: Derivative using FFT properties

MATLAB: Manual spectrogram creation without using spectrogram command

Related Question