【MATLAB】周波数が時間変化する信号を考える【瞬時周波数】

# はじめに
ある信号の周波数の時間変化を調べたい欲求がある人がいると思うが，そういうのを調べるのを時間周波数解析と言う．
時間周波数解析の基本として，今回は周波数の時間変化について解説し，代表的な周波数変調信号である線形チャープ信号や，任意の数式の周波数変調信号を作成し，`matlab`で確認してみることにする．

# 周波数が一定の信号を考える
周波数が時間変化する信号を考える前に，**周波数が時間変化しない**信号を考える．

## 1秒間に1回転する信号
```math:
y = \sin ( x )\\(0 \le x < 2 \pi)
```


この式は学校で習ったと思われるが，$x=2 \pi$で$y$が$1$回転して元に戻る数式である．
これを以下のように書き換えてみる．


```math:
y = \sin ( 2 \pi t )\\(0 \le t < 1)
```
こう書くと今度は$t=1$になったら$1$回転することになる．
$t$の単位を秒とすると，上の式は$1$秒間で$1$回転し，これが$1[\rm{Hz}]$の信号である．

これは円周上を回転する点であり，図にすると下のようになる．
<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/0e76ccfb-f4ba-7e9b-5b71-adb1a935add1.gif">
<details><summary>GIF画像の実装は↓</summary><div>

```matlab:
clc
clear all
close all
n0 = 0:0.01:2*pi;
F = 1;
t = linspace(0,1,20);
nImages = length(t);
fig = figure;
% アニメーションを1コマずつ作成する
for idx = 1:nImages
    plot(sin(n0),cos(n0),'LineWidth',3);
    hold on
    x = [0 cos(2*pi*F*t(idx))];
    y = [0 sin(2*pi*F*t(idx))];
    plot(x,y,'LineWidth',3);
    hold off
    title(['t = ' num2str( t(idx),'%0.5f\n') '[s]' ]);
    drawnow
    axis([-1,1,-1,1]);
    axis square
    frame = getframe(fig);
    im{idx} = frame2im(frame);
end
close;


filename = [num2str(F) 'Hz.gif'];
for idx = 1:nImages-1
    [A,map] = rgb2ind(im{idx},256);
    if idx == 1
        imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.05);
    else
        imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.05);
    end
end

```

参考：
https://jp.mathworks.com/help/matlab/ref/imwrite.html

</div></details>

次に，

```math:
y = \sin ( 2 \pi 2 t )
```

と書けば，$y$は$1$秒間に$2$周することになり，$2[\rm{Hz}]$の信号となる．（下図参照）

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/f9a1c042-52d8-3c6b-adb5-f489dae244dd.gif">





## 周波数一定の信号
以上の話をまとめると，周波数$f_c[\rm{Hz}]$（一定）で振動する信号は，

```math:
y = \sin ( 2 \pi f_c t )
```

と書ける．





# 周波数が時間変化する信号を考える

**周波数が時間変化する**とは，上の考え方で行くと，**円周上を動く点の速度が速くなったり遅くなったりする**ということになる．
図にすると以下のような感じ．

![fm01png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/39e215b9-3654-c483-07be-a97cf710fec2.png)

これは周波数が$1[\rm{Hz}]$から$2[\rm{Hz}]$に時間変化した例である．

なお，周波数が変化することを**変調**と言い，周波数が線形に変調する信号を**線形チャープ信号**や，**LFM信号（Linear Frequency Modulation signal，線形周波数変調信号）**と言ったりする．

## 瞬時周波数
さて，周波数が時々刻々と変化する場合において，ある時刻における一瞬の周波数はどのようにして求めたらいいだろうか．
これは**位相を時間で微分する**と求められ，ある時刻の瞬間の周波数が**瞬時周波数**である．

### 周波数一定の信号の瞬時周波数
再度，周波数$f_c$（一定）の信号で考えてみる．

```math:
y = \sin ( 2 \pi f_c t )
```

位相$\phi (t)=2 \pi f_c t$を時間微分して，

```math:
f_i(t) = \frac{1}{2\pi}\frac{d}{dt}  \phi (t) =\frac{1}{2\pi}\frac{d}{dt} ( 2 \pi f_c t ) = f_c \;\;({\rm constant})
```
周波数一定の信号の瞬時周波数$f_i=f_c$（一定）となり，期待通りである．


### 線形チャープ信号の瞬時周波数

$t=0[\rm{s}]$で$1[\rm{Hz}]$の信号が，周波数が線形に増大していき，$t=10[\rm{s}]$で$10[\rm{Hz}]$になる線形チャープ信号を考える．

$t$と$f_i$のグラフを以下に示す．
<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/9ee5bb60-7d32-7b0a-d344-2befbf800aee.png">


この数式は直線のY切片と傾きから，

```math:
f_i(t)= 0.9t+1
```

とすぐに求められる．


## 瞬時周波数から位相を求める
位相を時間で微分すると瞬時周波数が求められたので，逆に考えると，**瞬時周波数を時間で積分すると位相が求められる**ことになる．

```math:
\begin{align}
\phi (t)&= 2 \pi \int f_i dt\\
&=2 \pi \int (0.9t+1) dt\\
&=2 \pi  \left( \frac{0.9}{2}t^2+t \right) + \phi _0 \\
\end{align}
```

積分定数の$\phi _0$は初期位相である．

これを$\sin()$の中に入れて，

```math:
\begin{align}
y =\sin \left(2 \pi  \left( \frac{0.9}{2}t^2+t \right) + \phi _0 \right) \\
\end{align}
```

とすれば所望の線形チャープ信号が完成である．

### matlabで確認

```matlab:
figure
Fs = 1000;
t=0:1/Fs:10;
y=sin(2*pi*(0.9/2*t.^2 + 1*t));%初期位相は0とした
plot(t,y);
ylim([-1 1])
ylabel('振幅')
xlabel('t[s]')
```

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/60873986-b241-5230-1944-a3d9fdabae17.png">
という感じになり，sin波が時間経過でギューっと圧縮されており，確かに周波数がだんだんと高くなっているように見える．


### スペクトログラム
基本的な時間周波数解析として，短時間フーリエ変換を使って信号のスペクトログラムを求めてみる．

```matlab:
Nw = 2048;
Nol = round(Nw/2);
NFFT = Nw;
spectrogram(y, Nw, Nol, NFFT, Fs, 'yaxis');
ylim([0, 10])
```

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/cc4cbb24-a889-4e91-6724-cb137382c57c.png">

元の数式のグラフ（再掲）↓と比較して，同じような結果が得られた．

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/9ee5bb60-7d32-7b0a-d344-2befbf800aee.png">

## chirp関数
線形チャープ信号を手実装してみたが，`matlab`には`chirp`というチャープ信号を生成する関数が標準で備わっている．

[スイープ周波数の余弦信号 - MATLAB chirp - MathWorks 日本](https://jp.mathworks.com/help/signal/ref/chirp.html)

`chirp`は$\cos$で実装されているので，上の式を書き換えて以下のようにする．

```math:
\begin{align}
y =\cos \left(2 \pi  \left(\frac{0.9}{2}t^2+t\right) + \phi _0 \right) \\
\end{align}
```

```matlab:
figure
Fs = 1000;
t=0:1/Fs:10;

y=cos(2*pi*(0.9/2*t.^2 + 1*t));% 手実装

y2 = chirp(t,1,10,10);% matlab標準

plot(t,y);
hold on
plot(t,y2);
ylim([-1 1])
ylabel('振幅')
xlabel('t[s]')
legend('手実装', 'chirp関数')

% 確認
fprintf('norm=%f\n', norm(y-y2));% norm=0となれば両者は一致する
```
結果は下のようになり，ノルムは0となり，手実装したものと`matlab`標準の結果が一致した．$\sin$で作ったものと見比べると位相が$90$度だけずれていることも分かる．
<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/9076b116-30dd-2e53-f46f-1d0834871936.png">


## 任意のチャープ信号を作ってみる

大体のチャープ信号は`chirp`を使えば作れちゃうので，今度は`chirp`では作れない瞬時周波数を持つ信号を作成してみる．
瞬時周波数を以下のように定義する．

```math:
f_i(t) = 500 \sin (2 \pi t) + 500 t + 3000
```
瞬時周波数のグラフは以下となる．

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/e1c23e83-05e7-556a-5923-0d6e2bbe04b2.png">

時間$t=0[{\rm s}]$で$f=3000[{\rm Hz}]$，そこからほわんほわんと周波数が高くなっていく信号である．

先ほどと同様に，$f_i$を積分して位相を求める．

```math:
\begin{align}
\phi (t)&= 2 \pi \int f_i dt\\
&=2 \pi \int (500 \sin (2 \pi t) + 500 t + 3000) dt\\
&=2 \pi  \left( -\frac{500}{2 \pi} \cos (2 \pi t) + 250 t^2 + 3000t \right) + \phi _0\\
\end{align}
```

これを$\sin()$の中に入れて，

```math:
\begin{align}
y =\sin \left( 2 \pi  \left( -\frac{500}{2 \pi} \cos (2 \pi t) + 250 t^2 + 3000t \right) + \phi _0  \right)\\
\end{align}
```

時間波形は下のようになって，これだけだときしめんになっててよく分からない．
<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/83ceb87b-3f60-81dc-5531-45d61a350c8b.png">

2021年11月20日追記：上のきしめんだと不親切だと思ったので，ざっくりとした図を描いてみた（下図）．
時間周波数グラフと時間振幅グラフの対応がイメージ出来るようになっただろうか．
![ponti.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/e21574f5-5c27-4ce1-6f61-d2ccb35a1283.png)


さて，スペクトログラムを求めてみると，確かに期待した結果となった．
<img width="650" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/5e335e91-dd5a-5384-33b0-81df12b0721c.png">


#補足（物理の速度と瞬時周波数を比較する）
例題の瞬時周波数は以下であった．

```math:
f_i(t)= 0.9t+1
```
これに$2 \pi$を掛けたものが**瞬時位相**である．

```math:
\phi _i(t)= 2 \pi (0.9t+1)
```

これを以下のように書き換える．

```math:
v(t)=a t + v_0
```

こう書くと，これは物理の授業で習った物体の速度である．
物体の移動距離は以下のように書ける．

```math:
x(t)=\int v(t)dt = \frac{1}{2} a t^2 + v_0 t
```

まとめると，以下の表になる．

|                  | 速度の世界 |    周波数の世界    |
|:-----------------|:----------------------------------:|:------------------:|
|            |     距離 $x(t) = \frac{1}{2} a t^2 + v_0 t$ |      位相 $\phi(t) = 2 \pi  \left( \frac{0.9}{2}t^2+t \right)$  |
| $1$回微分          |      速度$v(t)=a t + v_0$ |瞬時位相$f_i(t)= 2 \pi (0.9t+1)$      |
| $2$回微分          |      加速度$a$ |$\frac{d}{dt}f_i(t)= 2 \pi 0.9$      |

加速度に対応する，瞬時周波数（瞬時位相）を微分したものがどういう名前で定義されているかは分からなかった．


https://twitter.com/dokagui_tairan/status/1461876525423370240?s=20

以上