MATLAB Language
Преобразования Фурье и обратные преобразования Фурье

Синтаксис

1. Y = fft (X)% вычисляет БПФ вектора или матрицы X, используя длину преобразования по умолчанию 256 (будет подтверждено для версии)

2. Y = fft (X, n)% вычисляет БПФ X, используя n в качестве длины преобразования, n должно быть числом на 2 степени. Если длина X меньше n, то Matlab автоматически вставит X с нулями, так что длина (X) = n

3. Y = fft (X, n, dim)% вычисляет БПФ X, используя n в качестве длины преобразования по размеру dim (может быть 1 или 2 для горизонтальной или вертикальной соответственно)

4. Y = fft2 (X)% Вычислить 2D БПФ X

5. Y = fftn (X, dim). Вычислить dim-мерный БПФ X относительно вектора размерностей dim.

6. y = ifft (X)% вычисляет Обратный БПФ X (который является матрицей / вектором чисел) с использованием 256 Transform Length

7. y = ifft (X, n)% вычисляет IFFT X, используя n в качестве длины преобразования

8. y = ifft (X, n, dim)% вычисляет IFFT X, используя n в качестве длины преобразования по размеру dim (может быть 1 или 2 для горизонтальной или вертикальной соответственно)

9. y = ifft (X, n, dim, 'симметричный).% Симметричный параметр заставляет ifft трактовать X как сопряженное симметричное по активной размерности. Эта опция полезна, когда X не является точно сопряженной симметричной, просто из-за ошибки округления.

10. y = ifft2 (X)% Вычислить обратный 2D ft X

11. y = ifftn (X, dim). Вычислить обратное dim-мерное fft пространства X.

параметры

замечания

ifft(fft(X)) = X

Внедрение простого преобразования Фурье в Matlab

Преобразование Фурье, вероятно, является первым уроком в цифровой обработке сигналов, это приложение повсеместно, и это мощный инструмент, когда дело доходит до анализа данных (во всех секторах) или сигналов. Matlab имеет набор мощных наборов инструментов для преобразования Фурье. В этом примере мы будем использовать преобразование Фурье для анализа основного синусоидального сигнала и генерировать то, что иногда называют периодограммой с использованием БПФ:

%Signal Generation
A1=10;                % Amplitude 1
A2=10;                % Amplitude 2
w1=2*pi*0.2;          % Angular frequency 1
w2=2*pi*0.225;        % Angular frequency 2
Ts=1;                 % Sampling time
N=64;                 % Number of process samples to be generated
K=5;                  % Number of independent process realizations
sgm=1;                % Standard deviation of the noise
n=repmat([0:N-1].',1,K);             % Generate resolution
phi1=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 1
phi2=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 2
x=A1*sin(w1*n*Ts+phi1)+A2*sin(w2*n*Ts+phi2)+sgm*randn(N,K);   % Resulting Signal

NFFT=256;            % FFT length
F=fft(x,NFFT);       % Fast Fourier Transform Result
Z=1/N*abs(F).^2;     % Convert FFT result into a Periodogram

Обратите внимание, что дискретное преобразование Фурье реализовано быстрым преобразованием Фурье (fft) в Matlab, оба будут давать тот же результат, но FFT - это быстрая реализация DFT.

figure
w=linspace(0,2,NFFT);
plot(w,10*log10(Z)),grid;
xlabel('w [\pi rad/s]')
ylabel('Z(f) [dB]')
title('Frequency Range: [ 0 ,  \omega_s ]')

Обратные преобразования Фурье

Одним из основных преимуществ преобразования Фурье является его способность обратно обратно в Time Domain без потери информации. Рассмотрим тот же Сигнал, который мы использовали в предыдущем примере:

A1=10;                % Amplitude 1
A2=10;                % Amplitude 2
w1=2*pi*0.2;          % Angular frequency 1
w2=2*pi*0.225;        % Angular frequency 2
Ts=1;                 % Sampling time
N=64;                 % Number of process samples to be generated
K=1;                  % Number of independent process realizations
sgm=1;                % Standard deviation of the noise
n=repmat([0:N-1].',1,K);             % Generate resolution
phi1=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 1
phi2=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 2
x=A1*sin(w1*n*Ts+phi1)+A2*sin(w2*n*Ts+phi2)+sgm*randn(N,K);   % Resulting Signal

NFFT=256;            % FFT length
F=fft(x,NFFT);       % FFT result of time domain signal

Если мы откроем F в Matlab, мы найдем, что это матрица комплексных чисел, вещественная часть и мнимая часть. По определению, чтобы восстановить исходный сигнал временной области, нам нужны как Real (который представляет вариацию Magnitude), так и Imaginary (который представляет изменение фазы), поэтому, чтобы вернуться во временную область, можно просто хотеть:

TD = ifft(F,NFFT);   %Returns the Inverse of F in Time Domain

Обратите внимание, что возвращаемый TD будет длиной 256, поскольку мы устанавливаем NFFT на 256, однако длина x равна только 64, поэтому Matlab будет заполнять нули до конца TD-преобразования. Так, например, если NFFT равнялся 1024, а длина равнялась 64, то возвращаемое TD было равно 64 + 960 нулей. Также обратите внимание, что из-за округления с плавающей запятой вы можете получить что-то вроде 3.1 * 10e-20, но для общего назначения: для любого X ifft (fft (X)) равно X с точностью до округлой ошибки.

Скажем на мгновение, что после преобразования мы сделали что-то и остались только с РЕАЛЬНОЙ частью БПФ:

R = real(F);         %Give the Real Part of the FFT
TDR = ifft(R,NFFT);  %Give the Time Domain of the Real Part of the FFT

Это означает, что мы теряем мнимую часть нашего БПФ, и поэтому мы теряем информацию в этом обратном процессе. Чтобы сохранить оригинал без потери информации, вы всегда должны сохранять мнимую часть БПФ с помощью imag и применять свои функции как к обеим, так и к реальной части.

figure
subplot(3,1,1)
plot(x);xlabel('time samples');ylabel('magnitude');title('Original Time Domain Signal')
subplot(3,1,2)
plot(TD(1:64));xlabel('time samples');ylabel('magnitude');title('Inverse Fourier Transformed - Time Domain Signal')
subplot(3,1,3)
plot(TDR(1:64));xlabel('time samples');ylabel('magnitude');title('Real part of IFFT transformed Time Domain Signal')

Изображения и многомерные ФТ

В медицинской визуализации, спектроскопии, обработке изображений, криптографии и других областях науки и техники часто бывает так, что хочется вычислить многомерные преобразования Фурье изображений. Это довольно просто в Matlab: (многомерные) образы - это просто n-мерные матрицы, в конце концов, и преобразования Фурье являются линейными операторами: одно просто итеративно преобразование Фурье по другим измерениям. Matlab обеспечивает fft2 и ifft2 чтобы сделать это в 2-d, или fftn в n-измерениях.

Одна потенциальная ошибка заключается в том, что преобразование Фурье изображений обычно показано «ориентированным по центру», т. Е. С началом k-пространства в середине изображения. Matlab предоставляет команду fftshift для правильной замены местоположения DC-компонентов преобразования Фурье. Это упорядочивающее обозначение значительно упрощает выполнение общих методов обработки изображений, один из которых проиллюстрирован ниже.

Нулевое заполнение

Один «быстрый и грязный» способ интерполировать небольшое изображение на больший размер - преобразовать его в Фурье, наложить преобразование Фурье на нули и затем принять обратное преобразование. Это эффективно интерполирует между каждым пикселем с основанной на сердцевиной базовой функцией и обычно используется для масштабирования медицинских изображений с низким разрешением. Начнем с загрузки примера встроенного изображения.

%Load example image 
I=imread('coins.png'); %Load example data -- coins.png is builtin to Matlab 
I=double(I); %Convert to double precision -- imread returns integers 
imageSize = size(I); % I is a 246 x 300 2D image

%Display it
imagesc(I); colormap gray; axis equal; 
%imagesc displays images scaled to maximum intensity

Теперь мы можем получить преобразование Фурье I. Чтобы проиллюстрировать, что делает fftshift , давайте сравним два метода:

% Fourier transform 
%Obtain the centric- and non-centric ordered Fourier transform of I
k=fftshift(fft2(fftshift(I))); 
kwrong=fft2(I); 

%Just for the sake of comparison, show the magnitude of both transforms: 
figure; subplot(2,1,1); 
imagesc(abs(k),[0 1e4]); colormap gray; axis equal; 
subplot(2,1,2); 
imagesc(abs(kwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear). 

Теперь мы получили 2D FT примерного изображения. Чтобы нуль заполнить его, мы хотим взять каждое k-пространство, напасть на края нулями, а затем взять обратное преобразование:

%Zero fill
kzf = zeros(imageSize .* 2); %Generate a  492x600 empty array to put the result in 
kzf(end/4:3*end/4-1,end/4:3*end/4-1) = k; %Put k in the middle 
kzfwrong = zeros(imageSize .* 2); %Generate a  492x600 empty array to put the result in 
kzfwrong(end/4:3*end/4-1,end/4:3*end/4-1) = kwrong; %Put k in the middle 

%Show the differences again 
%Just for the sake of comparison, show the magnitude of both transforms: 
figure; subplot(2,1,1); 
imagesc(abs(kzf),[0 1e4]); colormap gray; axis equal; 
subplot(2,1,2); 
imagesc(abs(kzfwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear). 

На данный момент результат довольно ничем не примечателен:

Как только мы возьмем обратные преобразования, мы можем видеть, что (правильно!) Данные с нулевым заполнением обеспечивают разумный метод интерполяции:

% Take the back transform and view 
Izf = fftshift(ifft2(ifftshift(kzf))); 
Izfwrong = ifft2(kzfwrong);

figure; subplot(1,3,1); 
imagesc(abs(Izf)); colormap gray; axis equal; 
title('Zero-filled image'); 
subplot(1,3,2); 
imagesc(abs(Izfwrong)); colormap gray; axis equal;
title('Incorrectly zero-filled image'); 
subplot(1,3,3); 
imagesc(I); colormap gray; axis equal; 
title('Original image'); 
set(gcf,'color','w'); 

Обратите внимание, что размер заполненного нулем изображения вдвое больше, чем у оригинала. В каждом измерении может быть заполнено нулевое значение более чем в два раза, хотя, очевидно, это не приводит к произвольному увеличению размера изображения.

Советы, подсказки, 3D и дальше

Вышеприведенный пример выполняется для трехмерных изображений (как это часто генерируется методами медицинской визуализации или конфокальной микроскопии, например), но требует, fft2 замены fftn(I, 3) на fftn(I, 3) . Из-за несколько громоздкой природы написания fftshift(fft(fftshift(... несколько раз, довольно часто для определения функций, таких как fft2c локально, для обеспечения более легкого синтаксиса локально - например:

function y = fft2c(x)

y = fftshift(fft2(fftshift(x)));

Обратите внимание, что FFT является быстрым, но большие многомерные преобразования Фурье по-прежнему занимают время на современном компьютере. Кроме того, он по своей сути является сложным: величина k-пространства была показана выше, но фаза абсолютно необходима; переводы в области изображения эквивалентны фазовой рампе в области Фурье. Есть несколько более сложных операций, которые, возможно, пожелают сделать в области Фурье, такие как фильтрация высоких или низких пространственных частот (путем умножения их на фильтр) или маскирование дискретных точек, соответствующих шуму. Соответственно, имеется большое количество кода, созданного сообществом для обработки общих операций Фурье, доступных на основном сайте репозитория сообщества Matlab, File Exchange .