・追記
一応、念のためさらに調査を行ったところ下記の記事が死ぬほどわかりやすく、はっきりと説明してくれていました。英語の読める方は、ここをみたら安心してスペクトルの計算ができると、思いまーす！

なお、本文にはいろいろと勘違いがあるので、最後の方にまとめて追記しています。

（追記ここまで）

いい記事ではないですが、よく分からないことがあったのでメモ程度に。数式が得意な人は、下記の記事の方が分かりやすいと思います。式がいっぱい載ってます。私は数式はダメなので、この記事ではうすぼんやりとメモ程度の説明をします。

Pythonのfftpack.fftを利用して、パワースペクトルだか、パワースペクトル密度だか、どっちかを計算したかったのですが、そもそもfftの戻り値が何なのかすら知らなかったので、いろいろと調べてみました。

どうやら、fftpack.fftはWikipediaにある通り、離散フーリエ変換を行った結果をそのまま返します。

要するに上の式で計算した値です。この辺、これが普通なのかどうか私は知りません。というのも、Wikipediaの全く同じページに次のような定義があるからです。

ぱっと見違う式に見えますが、

ということなので、違いはNで割っているかどうかだけです。離散フーリエ変換といった場合、どちらが正式な定義なのか、それとも正式な定義などないのか？このあたりについて、知っている人がいたら教えてくれると助かります。

で、本題なのですが、Nで割らない式を使うと、周波数領域で見たときのY軸、信号の振幅だか、強度の値だかが信号の長さ、及びサンプリング・レートによって違ってきます。長いデータであれば、大きな値になるわけです（これは、数学音痴の私でも式を見たらすぐにわかります）。ここで、はたとわからなくなってしまいました。理解がうすぼんやりなので、データ数Nで割るという操作が果たして「正しい」のかと考えると、そうでもない気がします。

色々調べて回った結果、スペクトルの表現には二種類あり、一つはパワー・スペクトルと呼ばれるもので、こちらは時間に依存しないらしいです。一方で、エネルギー・スペクトルというのもあり、こちらは信号のエネルギーに比例するわけですから、長い信号であればでかくなるはずです。一方で、離散フーリエ変換を行う場合、時間とサンプリング・レートには（信号のデータ点数N）＝（信号の持続時間T）×（サンプリング・レート）という関係がありますが、エネルギー・スペクトルは飽くまでエネルギーなので、サンプリング・レートとは何の関係もないはずです。

以上の考察から、もし、信号の長さに依存しない、Nで割った方の定義を「離散フーリエ変換」とするならば、パワースペクトルはfftpack.fftの戻り値（フーリエ変換の係数）をNで割ったものにすればよく、エネルギー・スペクトルにするにはサンプリング・レートで割ればよい（データ数Nで正規化した後、時間Tを掛けているのと同じ）ということで、たぶん、これが正解だと思います。というわけで、今作っているライブラリの実装はそのようにしました。

一方で、なんちゃらスペクトルには、密度とつくものとつかないものがあります。この違いもよく分からないのですが、どうやら、密度とつくものは、確率過程と関連付けられるようです。

ってことは、何らかのデジタル信号の周波数領域における強度を計算したいだけなら、別に密度じゃなくてもいいのかなーと思っていたら、冒頭でリンクしたQuiitaの記事に

「元の時間信号に含まれる単位周波数 (1 Hz幅) 当たりのパワー (パワー = 単位時間あたりのエネルギー) を周波数の関数として表現したもの」(ONOSOKKI3 pp. 2–3)，つまりパワースペクトルP(k)P(k)を単位周波数で規格化したスペクトルであり、周波数分解能をΔfとすると...

とサラっと書かれていたので、そうでもないようです。当該の記事で紹介されている数式は次の通り：

PSD(k)というのがパワースペクトル”密度”で、Δfが周波数分解能（サンプリング・レート）なので、パワースペクトルを「サンプリングレートで正規化したもの」がパワースペクトル密度のようです。要するに、フーリエ変換の結果そのまんまのやつを、時間で割ったものがパワースペクトル、それをさらにサンプリングレートで正規化したものが、パワースペクトル密度というのがたぶん正解のようですので、コードで書くと次のようになります。

・追記
これは間違いでした。サンプル数（時間×サンプリング・レート）をさらにサンプリング・レートでもっかい割ったのが、パワースペクトル密度だそうです。もしかすると分野によるのかもしれませんが、DSPなどの分野では単位をきちんとする関係からかそんな風に説明されており、明確でわかりやすかったので、それを信用することにしました。

（追記ここまで）

import numpy as np
from scipy.fft import fft


T=10 # 信号の持続時間
fs=100 # サンプリングレート
dt=1./fs # 一サンプル当たりの時間経過
ts=np.arange(0,T,dt)
f=1.0 # サイン波の周期
sig=np.sin(2*np.pi*f*ts) # sine信号の生成
fft_sig=fft(sig) # フーリエ変換
ps=fft_sig/T # パワースペクトルを計算 （←間違い）
psd=ps/fs # パワースペクトル密度を計算　（←これは密度ではなく、ただのパワースペクトル）
psd_good=psd/fs # これが本当の「パワースペクトル密度」らしい。後日追記。

>>> np.max(np.abs(ps))
50.0
>>> np.max(np.abs(psd)) # 普通、これが欲しいと思うの（けれどもよくよく調べてみるとこっちではなかった）
0.5

というわけで、パワースペクトルは周波数分解能に依存します。この例ではdtが有限なので、これでいいんじゃないかなー？ただ、これだと、周波数解像度をいくらでもあげれるのなら、パワースペクトルの値もいくらでも大きくなりそうな気がしますが…まぁ、数式だと、このdtがlimとってゼロの極限として積分するので、そうはならないからOKなのだと思います。というわけで、コンピューターで離散的な計算をする場合、こうやってパワースペクトル（密度じゃないほう）を計算するのは間違いなのかもしれません（なぜなら、コンピューターで計算する限り、飽くまで有限のΔtやΔfでなければ計算が成り立たないから）。よって、電子計算機でパワースペクトルを計算する、という場合、最後のpsdに相当するもの（パワースペクトル密度もどき？）を計算しようとしているのが正しいのではないかというのがここでの結論です。

・追記
他の箇所の追記にある通りこの結論は間違いです。ここで説明しているのはパワースペクトルの計算方法であり、密度を計算する場合は更にサンプリング・レートで割るそうです。これが、世間一般でいう「パワースペクトル密度」の意味かどうかは知りませんが、色々と納得がいったので、ここではそういうこととして理解しました。なんだか、不思議の国のアリスを思い出すなぁ。

（追記ここまで）

ちなみに、エネルギー・スペクトル”密度”が、psd*Tで良いのかどうか？はちょっとわかりませんでした。理屈から言うとそうなるのではないかと思うのですが、このあたりも知っている人がいたら、教えていただけると助かります。

ちなみに、このコードだと、一部、数学的な定義とつじつまが合わなくなるところがある気がします。というのも、ホワイトノイズ（ここではシンプルに、正規分布を用います）の分散は、実はパワースペクトルと関係があり、パワースペクトルの二乗はホワイトノイズの分散に等しくなるはずだからです。

というわけで、一応、適当なノイズを生成して計算してみると、ちゃんとPSDの絶対値の二乗を足し合わせたものと、元系列の分散は等しくなるので最初はOKかと思ったのですが、実は、ホワイトノイズを離散データとして生成する場合、分散はサンプリングレートに比例しますので、それも加えて考えると、パワースペクトルを計算する場合は、上記のコードでpsdとしているものを、さらにサンプリングレートで割らなければつじつまが合わないとも考えられます。（後日判明しましたが、あわないです。合わないのでサンプリング・レートでさらに割るのが正しいようです）。

うーん。どうなんでしょう。データ数Nで割ったものを、さらにサンプリングレートで割るというのは、つじつま合わせにしてもおかしい気がしたので、とりあえず、今は深く考えないことにします。実際、この場合「元系列の分散」をどう解釈するのかが問題です（元系列というか、ナイキスト周波数でカットする前の高周波成分まで含んだ真のノイズ源の分散）。とりあえず、データから計算される分散についてはつじつまが合っているのですから、これでいいっちゃいい気もします（ホワイトノイズの分散とサンプリングレートの関係は、後で補正することにして、そういうものと割り切る）。

ただし、ホワイトノイズのエネルギー（スペクトルを積分だか、合計だかしたもの）はサンプリングレートによって変わります。なぜなら、ホワイトノイズは周波数全域にまたがっているため、ナイキスト周波数でカットすることで、見かけ上、エネルギーが小さくなるからです。そういうことまで考えると、やっぱり「ホワイトノイズのスペクトル」を計算するのなら、その辺の帳尻もあっていないとダメなのかなぁという気もします。と勝手に思っていたのですが、ホワイトノイズのサンプリングと言うのは実は結構分かりにくい問題で、Stack Exchangeなどでもかなりのトピックを見つけることができます。

https://dsp.stackexchange.com/questions/56809/question-about-sampling-white-noise/56813#56813
https://dsp.stackexchange.com/questions/55132/power-spectral-density-of-discrete-white-gaussian-noise-defined-with-variance-an

一つ目の記事によると、ある周波数fにおけるホワイトノイズのスペクトル密度はΔt→0の極限（サンプリング・レートがいくらでも高くできるとき）Δt2σ^2…って書いています…え？えぇ…？。ただ、ホワイトノイズの自己相関は、デルタ関数ですから（二つ目の記事）フーリエ変換すると、一定値を取ったものがf=―∞から+∞までずっと続くことになります。こう考えると、サンプリング・レートが大きくなると、当然、ナイキスト周波数も高くなるので、スペクトル密度の面積が広がるから、そうはならないような気がします。まぁ、こんなことは専門家にしか分からなそうですが、なんか間違ってたら教えて欲しす。（追記：これも後々考えてみると、それほどおかしくもないのか？なぜなら、fsでもっかいわるってそういうことなんだから。要するに、フーリエ変換の計算で使うbinに相当するΔfが小さくなるんだから、それに比例してbin一個当たりに含まれるスペクトルも小さくなる…よって、Δtが小さくなれば、周波数Δfあたりのホワイトノイズのスペクトル密度も小さくならないとおかしいよね？という意味なんですね、うすぼんやりと）

というわけで、ネット上をあちこち駆けずり回ったのですが、パワースペクトルとパワースペクトル「密度」の違いはうすぼんやりのままです。

パワースペクトルは雑音の強度を各周波数成分に分解して、単位周波数あたりの強度で表わしたものである。

とかって書いてある資料もあったりして、え？えぇ？ってなってしまいました。これって、パワースペクトル「密度」の定義のはずじゃあ…まだ、何か勘違いしているのだろうか。パワースペクトルと、パワースペクトル密度の違いが、またうすぼんやりになってしまいました。これは、分野によって違うのか、その日によって違うのか、この言葉を使った人の朝食によって変わるのか…。

・追記
あとあともう少し詳しく調べてみると、やっぱり上記の記事にあった小野測器版の定義？が正しそうということが分かりました。ただし、私の解釈は間違いで、冒頭で紹介した英文の記事によると、

Consider a continuous broadband signal that is sampled at different sample rates but with the same number of samples. At a lower sample rate, each bin represents a narrower frequency band so that the magnitude of the power spectrum will be lower. The PSD compensates for this.

だそうで、power spectral densityと言う語が意味するものは、フーリエ変換の結果（scipy.fft.fftの戻り値）を、データ数Nで割り算した後、さらにサンプリングレートで割ったものということになります。これで、ホワイトノイズの分散がおかしくなる問題にも決着がつき、すべて整合性が取れました！やったね。

（追記ここまで）

というわけで、調べられる限り調べた結果、私が理解できた限りでは、これがパワースペクトル、パワースペクトル密度の「コード表現」になります。あっている保証はないですが、参考がてらメモ書きとして残しておきます。

あー、そうそう、あとpythonのfftモジュールはscipy.fftが一番新しい、おすすめのやつらしいです。np.fftとかscipy.fftpackは時代遅れだそうで…。

というわけで、うすぼんやりな理解ですが、参考になればと思います。