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

11 次查看（过去 30 天）
non
non 2022-11-10
评论： non 2022-11-17
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;
  1 个评论
Naoya
Naoya 2022-11-14
コードを確認する限りでは、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 のベクトルの要素に相当します)

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Hernia Baby
Hernia Baby 2022-11-16
Ran in:
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
  4 个评论
显示 2更早的评论
Hernia Baby
Hernia Baby 2022-11-17
Ran in:
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
non 2022-11-17
Hernia Baby様
たびたびのご回答を誠にありがとうございます。自身では到底考えつけない方法で、大変参考にさせていただきました。また、具体的なプログラムをご教示いただき、作図まで行うことができました。
ありがとうございました。

请先登录，再进行评论。

更多回答（1 个）

non
non 2022-11-14
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方向(振幅値)はすべて異なる値をもつように考えており、計算が正しいのか判断できない状況です。
お手数をおかけします。よろしくお願いいたします。
  6 个评论
显示 4更早的评论
Naoya
Naoya 2022-11-16
matファイルの中身を確認してみましたが、変数PC1M の行列が 49次元 (214x407x1x1x...x1x49) となっておりましたので、 214x407x49 にまずは修正の必要があります。
これについては、 squeeze コマンドで不要な次元削減を行うことができます。
PC1M = squeeze(PC1M);
non
non 2022-11-17
Naoya 様
何度もご回答いただくこととなり、申し訳ありませんでした。教えていただいたsqueezeコマンドを使って、望んでいた解析ができました。本当にありがとうございました。

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 PCM 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!