z変換による周波数解析について

Question

>> g=zeros(1,1000);
>> g(1:3)=[1 1 1];
>> fx=fft(g);
>> plot(20.*log10(abs(fx(1:500))))

このように入力するとg=[1 1 1]のz変換が求められるのですが,なぜこのように

g=zeros(1,1000);
g(1:3)=[1 1 1];

とセットしないと正しい周波数分解能が得られないのですか?

あとこれがその時に出力されたfftのグラフになるのですが,たとえばg=[1 1]であればゼロ点がz=-1でω=πのとことで応答が0になり,fs/2のところにディップが来ると思うのですが,2次の場合,z=(-1+√3i)/2,(-1-√3i)/2となるのですが,このようなグラフになる意味が理解できません.(-1+√3i)/2のゼロ点は340Hzのディップを表しているのではないかと思っているのですが,そこからどのように考えれば(-1-√3i)/2のゼロ点の位置がこのようになることを理解できますか.
質問がわかりにくいと思うので,追加すると半径1の1の単位円で考えた時にg=[1 1]でz=-1の時は単位円上に存在するのでfs/2のところにディップが来て,2枚目のグラフのように横軸の500Hzで0になっているのですが,2次になるとz=(-1+√3i)/2,(-1-√3i)/2でなぜこのようなグラフになるのか理解できません.z=(-1+√3i)/2,(-1-√3i)/2の点は単位円上にないからディップが下がりきらないのでしょうか.

また左が元の入力信号で右が得られたg=[1 1 1]のz変換を元の入力信号と畳み込んで,fftしたグラフなのです.
g=[1 1 1]のz変換はLPFのように感じるのですが,低域が逆に下がっています.なぜかよくわかりません.
よろしくお願いいたします.

Comments

[ozwk]

z変換と言っていますが、コードは離散フーリエ変換です。
(まあz=exp(iω)としたz変換ではありますが。)

ということは最初の3サンプルが1で続く997サンプルが0の方形波信号をただ離散フーリエ変換して強度をプロットしているだけなのですがそのあたり大丈夫ですか？

[kaeruuuun]

コメントありがとうございます.
大丈夫です.

[ozwk]

> なぜこのように
> とセットしないと正しい周波数分解能が得られないのですか?

むしろどうやろうとしていたのですか？
g=zeros(1,1000);を抜かすとかですか？

[ozwk]

> また左が元の入力信号で右が得られたg=[1 1 1]のz変換を元の入力信号と畳み込んで,fftしたグラフなのです.

畳み込みについてもコード貼ってください

[kaeruuuun]

g=zeros(1,1000);を抜かしたり,配列の要素が合わないので実行はできないのですが,g(1:3)=[1 1 1];の3を2にしたままで実行した場合です.

[kaeruuuun]

畳み込みはconvを使いました.コードは以下に示すようにwavを読み込んで,wavと畳み込みしました.
[x fs]=audioread('12.wav');
y=conv(g,x);

Answer 1

g=zeros(1,1000);

g(1:3)=[1 1 1];
とセットしないと正しい周波数分解能が得られないのですか?

g=zeros(1,1000);を抜かしたり,配列の要素が合わないので実行はできないのですが,g(1:3)=[1 1 1];の3を2にしたままで実行した場合です.

要素がN個のデータを離散フーリエ変換したらN個のデータが得られるので
gの要素がN=3個なら当然変換結果も3つです。

配列の要素が合わないのは
g(1:2)、つまりgの要素1から2までに1 1 1を代入しようとしています。
何がしたいのでしょうか。

2次になるとz=(-1+√3i)/2,(-1-√3i)/2でなぜこのようなグラフになるのか理解できません.z=(-1+√3i)/2,(-1-√3i)/2の点は単位円上にないからディップが下がりきらないのでしょうか.

(-1±√3i)/2は単位円上の点です。
2π/3と4π/3ですね。要するに3分割。

データ数が1000個なので333.333...(と666.666...)の位置です。

Comments

[kaeruuuun]

>>g(1:2)、つまりgの要素1から2までに[1 1 1](要素3つ)を代入しようとしています。
何がしたいのでしょうか。
とのことですが,たとえばg=zeros(1:1000)を抜かしたり,gの配列をうまく設定できなかった場合,周波数分解能がどうなるのか知りたいということです.
>>(-1±√3i)/2は単位円上の点です。
2π/3と4π/3ですね。要するに3分割。
これは2π/3が333.333の位置のディップで画像では表示されていないですが4π/3が666.666の位置で同じようにディップが存在するのでしょうか.
単位円で考えると0からπで一つ目の2π/3のディップがあり,πから-πにかけてまた山が上がって4π/3でディツプが発生するのでしょうか.

[jbpb0]

質問者さん

> g=zeros(1:1000)を抜かしたり,gの配列をうまく設定できなかった場合,周波数分解能がどうなるのか

上記の「周波数分解能」が、FFT結果の周波数の間隔の意味ならば、それは元データの全長さの逆数です

たとえば、1秒間隔に取得したデータが10個ある場合、周波数の間隔は1/10=0.1Hzなので、FFT結果の周波数は、順番に下記の通りになります
0, 0.1, 0.2, 0.3, 0.4, 0.5, -0.4, -0.3, -0.2, -0.1 Hz

1秒間隔で100個の場合は、周波数の間隔は同様に1/100=0.01Hzなので、FFT結果の周波数は、順番に下記の通り (長いので途中は省略)
0, 0.01, 0.02...0.48, 0.49, 0.5, -0.49, -0.48...-0.02, -0.01 Hz

元データが、最初の2～3個だけが値を持ち、それ以外は0の場合は、0のデータの数を変えて、データ総数が10個、100個、1000個となるデータを作り、それぞれのFFT結果を同じグラフに重ね書きしたら、FFT結果の周波数間隔が変わることをちゃんと考慮してグラフを書けば、グラフは重なります

元々の質問文をちゃんと読んでないので、質問者さんの疑問とずれたコメントになってるかもしれません
もし外したコメントなら、スルーしてください

[kaeruuuun]

コメントありがとうございます.
>>元データが、最初の2～3個だけが値を持ち、それ以外は0の場合は、0のデータの数を変えて、データ総数が10個、100個、1000個となるデータを作り、それぞれのFFT結果を同じグラフに重ね書きしたら、FFT結果の周波数間隔が変わることをちゃんと考慮してグラフを書けば、グラフは重なります

たとえば元データが1000個で最初の2個だけ値を持って残りが0の場合と最初の3個だけ値を持って残りが0の場合では周波数分解能は元データの数が一緒なら変わらないですよね?

[jbpb0]

> 周波数分解能は元データの数が一緒なら変わらないですよね?

元データのサンプリングレート(私が書いた例では1秒に1回)も同じなら、そうです

[kaeruuuun]

わかりました.
ありがとうございます.