3次元配列をフーリエ​​変換し、座標ごとの​振​幅を取り出すには​どう​すればよいです​か?

data=214(x)×407(y)×49(t)の3次元配列を、xyの座標位置ごとにフーリエ変換したいです(214×407回分、別々にフーリエ変換したいです)。

また、計算される座標ごとの振幅の大きさを、各座標の濃淡で表現したいです。

最終的に、周波数ごとに振幅の大きさが分かるような画像を作りたいと考えています。なお、もとの49次元の画像には値のない画素(Nan)があります。

下記のコードを実行すると、フーリエ変換した変数Yが49のz軸の方向に、同じ値で計算されました(Y(:,:,1)とY(:,:,2)が同じ行列になりました。z軸の方向にフーリエ変換した結果を振幅で取り出したいのですが、その方法をご教示いただけないでしょうか。

前回こちらでご回答いただき、途中まで解決したものの、上記のことが分からず、重複のご質問をしております。よろしくお願いいたします。

L=49;
Fs=1;
n=2^nextpow2(L);
Y=fft(data,n,3);
P2 = abs(Y/L);
P1 = P2(:,:,1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;

[Naoya]

コードを確認する限りでは、xyの座標位置ごとに1次元のフーリエ変換をz方向に掛けていることはご所望の通り実施できている模様です。

P1(2:end-1) = 2*P1(2:end-1);

の部分は、通常、フーリエ変換を掛ける次元方向に適用するものになりますので、

P1(:,:,2:end-1) = 2*P1(:,:,2:end-1);

が所望になると思います。

また、

n=2^nextpow2(L);
Y=fft(data,n,3);

で 64点で fft を行っていますので、 Y の 3次元めのサイズは 64となりますので、

f は、

f = Fs*(0:(n/2))/n;

になります。

周波数ごとにxy座標の振幅値を求めるのであれば、 P(:,:,N) で求められます。 (ここで N は 周波数ビンとなり、 1 ～ 33 となり、それぞれ f のベクトルの要素に相当します)

Answer 1

[Hernia Baby]

5×5の cell データを作ります

clear,clc,close all;
Fs = 1000;
t  = (0:1/Fs:1-1/Fs)';
f  = linspace(10,490,25)

f = 1×25

10 30 50 70 90 110 130 150 170 190 210 230 250 270 290 310 330 350 370 390 410 430 450 470 490

for ii = 1:length(f)
    C{ii,1} = sin(2*pi*f(ii)*t) + 0.01.*randi([0 1],size(t));
end

cellの配置方法はあまり考えてません

reshapeで並べ替えます

C = reshape(C,[],sqrt(length(f)));

cellfunでFFTをし、振幅をとります

C_fft = cellfun(@(x) abs(fft(x)),C ,'UniformOutput' ,false);

対応する周波数を作成します

freq = (linspace(0,Fs,length(t)))';

今回は150Hz付近を抜き出します

f1 = 150;
idx = knnsearch(freq,f1);

各要素から120Hzの振幅をとります

C_map = cellfun(@(x) x(idx),C_fft);

マップを作成します

x = [1 width(C_map)];
y = [1 height(C_map)];

figure
hold on
imagesc(x,y,C_map)
axis([0 x(end)+1 0 y(end)+1])

title(sprintf('%iHz付近の各振幅',f1))
colorbar

[non]

Hernia Baby様

ご回答いただきまして、誠にありがとうございます。

ご教示いただきました方法をそのまま適用することができておらず、1点ご質問があります。

C(5×5cell)・C{1,1}(1000×1double)を自身のデータで作成したいです。

今持っているデータは、PC1M=214×407×49(49-D double)でして、これをPC1M(214×407cell)・PC1M{1,1}(49×1double)に変換することは可能でしょうか。PC1Mを一度分解してから値を取り出す方法ではうまくいかず、もし可能でしたら教えていただきたいです。

ご確認をよろしくお願いいたします。

[Hernia Baby]

可能です。

以下のようにすれば問題ないです。

% データ作成
Data = rand(214,407,49);

% サイズ確認
Dx = height(Data)

Dx = 214

Dy = width(Data)

Dy = 407

% セルに格納
for ii = 1:Dx
    for jj = 1:Dy
        C{ii,jj} = squeeze(Data(ii,jj,:));
    end
end

% サイズ確認
size(C)

ans = 1×2

214 407

% 中身のサイズ確認
size(C{1,1})

ans = 1×2

49 1

[Hernia Baby]

matデータがあるんですね

具体的には以下になります

load 'PC1M.mat'
Data = squeeze(PC1M);
size(Data)

ans = 1×3

214 407 49

Dx = height(Data);
Dy = width(Data);
% セルに格納
for ii = 1:Dx
    for jj = 1:Dy
        C{ii,jj} = squeeze(Data(ii,jj,:));
    end
end

% サイズ確認
size(C)

ans = 1×2

214 407

% 中身のサイズ確認
size(C{1,1})

ans = 1×2

49 1

[non]

Hernia Baby様

たびたびのご回答を誠にありがとうございます。自身では到底考えつけない方法で、大変参考にさせていただきました。また、具体的なプログラムをご教示いただき、作図まで行うことができました。

ありがとうございました。

Answer 2

[non]

Naoya様

ご回答いただき、誠にありがとうございます。

ご指摘いただいた方法に修正して、プログラムを動かしてみました。

まだ、計算される値が正しいのかについて確証が持てず、一つ、ご質問させていただけないでしょうか。

L=49;
Fs=1;
n=2^nextpow2(L);

Y=fft(data,n,3);
P2 = abs(Y/L);
P1 = P2(:,:,1:L/2+1);

P1(:,:,2:end-1) = 2*P1(:,:,2:end-1);
f = Fs*(0:(n/2))/n;


a= P1(:,:,1);
aa=P1(:,:,2);
aaa=P1(:,:,24);
aaaa=P1(:,:,25);

と計算した結果、a(214×407 double)とaaaa(214×407 double)およびaa(214×407 double)とaaa(214×407 double)が全く同じ行列となりました。また、aとaaaaの行列の2倍が、aaとaaaに等しくなりました。これは、どのように理解すればよいか、可能であれば教えていただけないでしょうか。

自身では、フーリエ変換後のz方向(振幅値)はすべて異なる値をもつように考えており、計算が正しいのか判断できない状況です。

お手数をおかけします。よろしくお願いいたします。

[Naoya]

さらに確認をさせていただきましたが、以下2点の修正は必要の模様です。

% nextpow2 は 2のべき乗の指数部を返すので以下が正しい記述となります

n= nextpow2(L) --> n=2^nextpow2(L);

% FFTの点数はL ではなくn点となっていますので以下がよろしいと思います

P1 = P2(:,:,1:L/2+1); --> P1 = P2(:,:,1:n/2+1);

[non]

Naoya様

ご回答いただき、誠にありがとうございます。

ご指摘いただいた箇所を変更して、現在以下のコードを考えております。

ただ、まだ25次元すべてで違う振幅値を計算できておらず、以下のようになる原因が分からない状況です。

---

a(214×407 double)とaaaa(214×407 double)およびaa(214×407 double)とaaa(214×407 double)が全く同じ行列となりました。また、aとaaaaの行列の2倍が、aaとaaaに等しくなりました。

---

もし、何かお気づきの点がありましたら、ご教示いただけないでしょうか。

初学者でご質問ばかりとなり、申し訳ありません。

ご確認をよろしくお願いいたします。

L=49;
Fs=1;
n=2^nextpow2(L);

Y=fft(PC1M,n,3);
P2 = abs(Y/L);
P1 = P2(:,:,1:n/2+1);

P1(:,:,2:end-1) = 2*P1(:,:,2:end-1);
f = Fs*(0:(n/2))/n;


a= P1(:,:,1);
aa=P1(:,:,2);
aaa=P1(:,:,24);
aaaa=P1(:,:,25);

[Naoya]

PC1M をrand関数で3次元行列を疑似的に作って確認してみましたが、 a,aa,aaa,aaaa はそれぞれ異なる 2次元行列となっておりました。

3次元行列 PC1M について、PC1M(:,:,1), PC1M(:,:,2), PC1M(:,:,3) ... にて何かお気づきの点などありますでしょうか？

[non]

Naoya様

たびたびのご回答を誠にありがとうございます。

PC1M(:,:,1), PC1M(:,:,2), PC1M(:,:,3) ...には、すべて同じデータが入っております。PC1Mの値が、49-D doubleとして表示されること(214×407×49 doubleでないこと)は何かエラーの原因になりうるでしょうか。

変換しているファイルを添付しております。

もし可能でしたら、ご確認をよろしくお願いいたします。

[Naoya]

matファイルの中身を確認してみましたが、変数PC1M の行列が 49次元 (214x407x1x1x...x1x49) となっておりましたので、 214x407x49 にまずは修正の必要があります。

これについては、 squeeze コマンドで不要な次元削減を行うことができます。

PC1M = squeeze(PC1M);

[non]

Naoya 様

何度もご回答いただくこととなり、申し訳ありませんでした。教えていただいたsqueezeコマンドを使って、望んでいた解析ができました。本当にありがとうございました。