MATLAB Language
フーリエ変換および逆フーリエ変換
サーチ…
構文
Y = fft(X)%は、ベクトルまたは行列XのFFTを、デフォルトのTransform Length 256(バージョンの確認が必要)を使用して計算します。
Y = fft(X、n)%は、Transform Lengthとしてnを使用するXのFFTを計算します.nは2乗に基づく数値でなければなりません。 Xの長さがnより小さい場合、MatlabはXに自動的に0を埋め込み、length(X)= n
Y = fft(X、n、dim)%は、次元dimに沿ってTransform Lengthとしてnを使用してXのFFTを計算します(水平または垂直それぞれ1または2になります)
Y = fft2(X)%Xの2D FFTを計算する
Y = fftn(X、dim)%次元のベクトルに対してdimの次元FFTを計算する。
y = ifft(X)%は、デフォルトの256のTransform Lengthを使って、XのFFTの逆数(これは行列の数/ベクトルです)を計算します
y = ifft(X、n)%は、Transform Lengthとしてnを使用してXのIFFTを計算します。
y = ifft(X、n、dim)%は、次元dimに対してTransform Lengthとしてnを使用してXのIFFTを計算します(それぞれ水平または垂直に対して1または2になります)
y = ifft(X、n、dim、 'symmetric')%Symmetricオプションは、Xをアクティブ次元に沿って共役対称として扱うifftを引き起こします。このオプションは、Xが正確に共役対称でない場合に便利です。丸め誤差のためだけです。
y = ifft2(X)%Xの逆2次元ftを計算する
y = ifftn(X、dim)%Xの逆dim次元fftを計算します。
パラメーター
パラメータ | 説明 |
---|---|
バツ | これはあなたの入力のTime-Domain信号です。数字のベクトルでなければなりません。 |
n | これはTransform Lengthと呼ばれるNFFTパラメータで、FFT結果の分解能と考えると、2の累乗(つまり64,128,256 ... 2 ^ N)の数値でなければなりません(MUST) |
暗い | これは、FFTを計算する次元です。水平方向にFFTを計算する場合は1を使用し、垂直方向にFFTを計算する場合は2を使用します。 - このパラメータは通常空白のままです。あなたのベクトルの方向を検出することができます。 |
備考
Matlab FFTは、大量のデータを処理することができる非常に並列化されたプロセスです。また、GPUを使って大きな利点を得ることもできます。
ifft(fft(X)) = X
上記の文は、丸め誤差が省略されている場合に真です。
Matlabで簡単なフーリエ変換を実装する
フーリエ変換はおそらくデジタル信号処理の最初の教訓です。アプリケーションはどこにでもあり、データ(すべてのセクタ)や信号を分析する際には強力なツールです。 Matlabには、フーリエ変換用の強力なツールボックスが用意されています。この例では、フーリエ変換を使用して基本正弦波信号を解析し、時々FFTを使用してピリオドグラムと呼ばれるものを生成します。
%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
Discrete Fourier TransformはMatlabの高速フーリエ変換(fft)によって実装され、どちらも同じ結果が得られますが、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 ]')
逆フーリエ変換
フーリエ変換の大きな利点の1つは、情報を失うことなく時間領域に逆戻りする能力です。前述の例で使用したのと同じ信号を考えてみましょう。
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
MatlabでF
を開くと、それは複素数の行列、実数部と虚数部であることがわかります。定義上、元のタイム・ドメイン信号を回復するには、Real(マグニチュード・バリエーションを表す)とImaginary(フェーズ・バリエーションを表す)の両方が必要なので、タイム・ドメインに戻るには、
TD = ifft(F,NFFT); %Returns the Inverse of F in Time Domain
私たちがNFFTを256に設定しているので返されるTDは長さ256になりますが、xの長さは64にすぎないので、MatlabはTD変換の最後に0を埋め込みます。たとえば、NFFTが1024で長さが64の場合、返されるTDは64 + 960のゼロになります。また、浮動小数点の丸めのために、3.1 * 10e-20のようなものが得られるかもしれないが、一般的な目的のために:任意のXに対して、ifft(fft(X))は丸め誤差内のXに等しい。
変換を行った後、私たちは何かをしただけで、FFTの本当の部分が残っていると言いましょう。
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
これは、FFTの虚数部分を失っていることを意味し、したがって、この逆のプロセスで情報を失いつつあります。情報を失うことなくオリジナルを保存するには、 imag
を使用してFFTの虚数部を保持し、関数を実部または実部のいずれかに適用する必要があります。
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')
画像と多次元FT
医用画像、分光、画像処理、暗号、および科学技術の他の分野では、画像の多次元フーリエ変換を計算することがしばしば望ましい。結局のところ、(多次元の)多次元画像はn次元の行列であり、フーリエ変換は線形演算子です。他の次元に沿ったフーリエ変換だけです。 MATLABが提供fft2
およびifft2
2次元でこれを行うために、又はfftn
N次元で。
1つの潜在的な落とし穴は、画像のフーリエ変換は、通常、「中心の順序付け」、すなわち画像の中央にk空間の原点が示されていることである。 Matlabはフーリエ変換のDC成分の位置を適切に入れ替えるためのfftshift
コマンドを提供します。この順序記法は、一般的な画像処理技術を実行することを実質的により容易にし、その1つを以下に示す。
ゼロ充填
小さい画像をより大きなサイズに補間する「速くて汚い」方法の1つは、フーリエ変換し、フーリエ変換をゼロで埋め込み、次いで逆変換を行うことである。これは、各画素間で、シンク形状の基底関数を効果的に補間し、一般に、低解像度医用画像データをアップスケールするために使用される。組み込みの画像の例を読み込んでみましょう
%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
fftshift
が何をするかを説明するために、2つの方法を比較しましょう:
% 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空間を取り、エッジを0で埋め込み、バックトランスフォームを取ります。
%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');
ゼロ充填画像サイズは元の画像サイズの2倍です。明らかにそうしてもイメージのサイズは任意に増やされるわけではありませんが、各次元で2の倍以上でゼロを塗りつぶすことができます。
ヒント、ヒント、3D、それ以降
上記の例では、(多くの場合、例えば、医療用画像化技術または共焦点顕微鏡によって生成された)3D画像にも当てはまるが、必要fft2
によって置換されるfftn(I, 3)
例えば、。 fftshift(fft(fftshift(...
何回か)、 fft2c
などの関数をローカルで簡単に定義できるようにするのはかなり一般的です。
function y = fft2c(x)
y = fftshift(fft2(fftshift(x)));
FFTは高速ですが、大規模で多次元のフーリエ変換は現代のコンピュータではまだ時間がかかることに注意してください。それはさらに本質的に複雑です.k空間の大きさは上に示されていますが、位相は絶対に重要です。画像領域における平行移動は、フーリエ領域における位相傾斜と等価である。フーリエ領域では、高空間周波数または低空間周波数をフィルタリングする(フィルタに掛ける)か、ノイズに対応する離散点をマスクするなど、はるかに複雑な操作がいくつかあります。これに対応して、Matlabの主要コミュニティリポジトリサイトであるFile Exchangeで利用可能な共通フーリエ演算を処理するためのコミュニティ生成コードが大量に存在します。