Here's some code to calculate and display the FFT. It's in Matlab, but you're using Octave, so you should be able to adapt it.
classdef RawSpectrumAnalyzer
% RawSpectrumAnalyzer - Computes average horizontal and vertical power spectra
% from a raw image file (e.g., .dng) using rawread.
properties
fileName = '_DSF0027.dng' % Path to the raw image file (.dng, .cr2, etc.)
planeNumber = 1 % Raw plane index to extract
rawImage % Raw image data
imageSize % [rows, cols] of the selected plane
blackLevel = 64;
upperYlim = -10;
end
methods
function obj = RawSpectrumAnalyzer(fileName, planeNumber)
% Constructor
if nargin > 0
obj.fileName = fileName;
end
if nargin > 1
obj.planeNumber = planeNumber;
end
obj = obj.loadRawImage();
end
function obj = loadRawImage(obj)
% Load raw image and extract the specified plane
raw = rawread(obj.fileName);
raw = raw - obj.blackLevel;
if ndims(raw) == 2
img = raw;
elseif ndims(raw) == 3
if obj.planeNumber > size(raw, 3)
error('planeNumber exceeds number of planes in raw image.');
end
img = raw
, :, obj.planeNumber);else
error('Unsupported raw image format.');
end
obj.rawImage = double(img);
obj.imageSize = size(img);
end
function [horizontalSpectrum, verticalSpectrum] = computeSpectra(obj)
% Compute average horizontal and vertical spectra
img = obj.rawImage;
[rows, cols] = size(img);
obj.imageSize(1) = rows;
obj.imageSize(2) = cols;
horizontalSpectrum = zeros(1, cols);
for r = 1:rows
F = fft(img(r,
);P = abs(F).^2;
horizontalSpectrum = horizontalSpectrum + P;
end
horizontalSpectrum = horizontalSpectrum / rows;
verticalSpectrum = zeros(1, rows);
for c = 1:cols
F = fft(img
, c));P = abs(F).^2;
verticalSpectrum = verticalSpectrum + P.';
end
verticalSpectrum = verticalSpectrum / cols;
end
function plotSpectra(obj)
[horz, vert] = obj.computeSpectra();
rows = obj.imageSize(1);
cols = obj.imageSize(2);
horz = horz(1:floor(cols/2));
vert = vert(1:floor(rows/2));
horz = log2(max(horz, 1e-12) / max(horz));
vert = log2(max(vert, 1e-12) / max(vert));
f_h = linspace(0, 0.5, numel(horz));
f_v = linspace(0, 0.5, numel(vert));
figure;
subplot(2,1,1);
plot(f_h, horz);
title('Average Horizontal Spectrum');
xlabel('Spatial Frequency (cycles/pixel)');
ylabel('Relative Power (stops)');
drawnow; % Force plot rendering to get Y limits
yl = ylim;
ylim([yl(1), obj.upperYlim]);
subplot(2,1,2);
plot(f_v, vert);
title('Average Vertical Spectrum');
xlabel('Spatial Frequency (cycles/pixel)');
ylabel('Relative Power (stops)');
drawnow;
yl = ylim;
ylim([yl(1), obj.upperYlim]);
end
end
end
