FFT結果を片側スペクトルを行い保存する方法

[x,fs]=audioread("gum70-25.wav");

f = 50:50:500;

Fs = 44100;

dt=1/Fs;

duration=(length(x)-1)*(1/fs);

t = (0:1/Fs:duration)';

n = floor(length(t)/Fs);

w = hamming(Fs);

for ii = 1:n

z{ii} = fft(w.*x(Fs*(ii-1)+1:Fs*ii));

z1{ii}=abs(z{ii})

z2{ii}= z1{ii}(10001:20000,:);

writematrix(z2{ii},sprintf('gum70-25_fftsave_%03d.csv',ii))

end

wavファイルを用いてワークスペースにベクトル化を行い、その後FFTを44100点ずつ行いました。FFTを行った後に振幅スペクトルを計算しその後1ok~20kHzの部分のみを抜き出そうとしたのですが、うまくできませんでした。

簡単なことではありますがご教授おねがいいたします。

[Hernia Baby]

うまくいかない部分は出力でしょうか？バンドパスの部分でしょうか？

[一秀 近藤]

返信ありがとうございます。

バンドパスの部分です。

振幅スペクトルの仕方がまちがっているのかなと考えています。信号の片側振幅スペクトルを計算し、その中の10k-20kのみをデータとして保存したいと考えてます。

Answer 1

[Hernia Baby]

stftを使って振幅の片側スペクトルを求めてます。

今回はサンプリング周波数10,000Hz, 窓の大きさを1,000とします。

clear,clc;
Fs       = 10000;
winsize  = 1000;
Bandpath = [180 220];

サンプルデータを作ります。

t = (0:1/Fs:3-1/Fs).';
x = 0.5*cos(2*pi*t*200) + 0.1*randn(size(t));

figure
plot(t,x,'Color',[.4 .4 .4])
xlabel('時間[sec')
ylabel('信号')

stftで各パラメタを設定します。

今回窓は箱型にしています。

[s,f,t] = stft(x,Fs,'OverlapLength',0,'FFTLength',winsize,'FrequencyRange','onesided','Window',rectwin(winsize));

ここでバンドパス(というか周波数の範囲を決めただけ)を行います。

振幅については正規化します。

idx = f >= Bandpath(1) & f <= Bandpath(2);
s_b = abs(s(idx,:))./(winsize/2);
f_b = f(idx);

可視化します。

figure
heatmap(t,f_b,s_b)
xlabel('時間[sec]')
ylabel('周波数[Hz]')
colorbar
colormap(jet(10))
clim([0 1])