功率譜和頻譜：
功率譜：信號自相關后FFT
頻譜：信號直接FFT
功率譜：
信號的傳播都是看不見的，但是它以波的形式存在着，這類信號會產生功率，單位頻帶的信號功率就被稱之為功率譜。它可以顯示在一定的區域中信號功率隨着頻率變化的分布情況。
功率譜可以從兩方面來定義：
一個 是自相關函數的傅立葉變換；（維納辛欽定理）
另一個 是時域信號傅氏變換模平方然后除以時間長度。（來自能量譜密度） 根據parseval定理，信號傅氏變換模平方被定義為能量譜，能量譜密度在時間上平均就得到了功率譜。 頻譜：
頻譜是常常指信號的Fourier變換。 （1-7 作者：Yorkxu）轉載的理解：
（1）信號通常分為兩類：能量信號和功率信號；
能量信號：又稱能量有限信號，是指在所有時間上總能量不為零且有限的信號。 功率信號：它的能量為無限大，它對通信系統的性能有很大影響，決定了無線系統中發射機的電壓和電磁場強度。
（2）一般來講，能量信號其傅氏變換收斂（即存在），而功率信號傅氏變換通常不收斂（當然，若信號存在周期性，可引入特殊數學函數（Delta）表征傅氏變換的這種非收斂性）；
（3）信號是信息的搭載工具，而信息與隨機性緊密相關，所以實際信號多為隨機信號，這類信號的特點是狀態隨機性隨時間無限延伸，其樣本能量無限。換句話說，隨機信號（樣本）大多屬於功率信號而非能量信號，它並不存在傅氏變換，亦即不存在頻譜；
（4）若撇開搭載信息的有用與否，隨機信號又稱隨機過程，很多噪聲屬於特殊的隨機過程，它們的某些統計特性具有平穩性，其均值和自相關函數具有平穩性。對於這樣的隨機過程，自相關函數蛻化為一維確定函數，前人證明該確定相關函數存在傅氏變換；
（5）能量信號頻譜通常既含有幅度也含有相位信息；幅度譜的平方（二次量綱）又叫能量譜（密度），它描述了信號能量的頻域分布；功率信號的功率譜（密度）描述了信號功率隨頻率的分布特點（密度：單位頻率上的功率），業已證明，平穩信號功率譜密度恰好是其自相關函數的傅氏變換。對於非平穩信號，其自相關函數的時間平均（對時間積分，隨時變性消失而再次退變成一維函數）與功率譜密度仍是傅氏變換對；
（6）實際中我們獲得的往往僅僅是信號的一段支撐，此時即使信號為功率信號，截斷之后其傅氏變換收斂，但此變換結果嚴格來講不屬於任何“譜”（進一步分析可知它是樣本真實頻譜的平滑：卷積譜）；
（7）對於（6）中所述變換若取其幅度平方，可作為平穩信號功率譜（密度）的近似，是為經典的“周期圖法”；
補充：（8） 一個信號的頻譜,只是這個信號從時域表示轉變為頻域表示,只是同一種信號的不同的表示方式而已；而功率譜是從能量的觀點對信號進行的研究,其實頻譜和功率譜的關系歸根揭底還是信號和功率,能量等之間的關系。
譜估計 功率譜估計一般分成兩大類： 經典譜估計，也稱為非參數譜估計。 現代譜估計，也稱為參數譜估計。
clc;close all;clear all;
[s, fs] = audioread('hello.wav');%s=166912*1;fs=44100
%命令說明：[y,Fs] = audioread(filename)：returns sampled data, y, and a sample rate for that data, Fs.
%sound(y,fs) % set analysis parameters, pre-emphasise and windowing %根據話音位置，取5000可以把話音主要部分加入其中
figure; subplot(2,1,1); plot(s); ylabel('振幅'); xlabel('Time (n)'); title('原信號的時域'); Xk=fft(s);
subplot(212);plot(abs(Xk)); xlabel('\omega/\pi');ylabel('e^j^\omega');title('原信號的頻域');
%===============================================================
N = 4000;
Nfft = 4000;
n0 = 10000;
x = s(n0 : n0+N-1);%從10000個點截到14000，共4000個點
x1 = filter([1 -0.97], 1,x); %預加重 濾波器 %filter：Y = FILTER(B,A,X) ，輸入X為濾波前序列，Y為濾波結果序列，B/A 提供濾波器系數，B為分子， A為分母
%整個濾波過程是通過下面差分方程實現的： %a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + … + b(nb+1)*x(n-nb)-a(2)*y(n-1) - a(3)*y(n-2) + … + a(nb+1)*y(n-nb)
%作用：消除6dB/oct（分貝/倍頻程）的跌落，使語音信號的頻譜變得平坦。
w = (window('hamming', N)); xw = x1 .* w; %加窗 % Estimate PSD of the short-time segment
Sxw = fft(xw, Nfft);
Sxdb = 20*log10(abs(Sxw(1 : Nfft/2))) - 10*log10(N); %Sxdb 功率譜：時域fft取模平方后除以信號的長度 轉換成db
figure;
subplot(311);
plot(x);ylabel('振幅'); xlabel('Time (n)'); title('截取信號的時域信號');
subplot(312);
plot(x1);ylabel('振幅');xlabel('Time(n)');title('截取信號通過filter預處理后的時域信號');
subplot(313);
plot(xw);ylabel('振幅'); xlabel('Time (n)'); title('截取信號加漢明窗后的時域信號');
%============================================================================
figure;
subplot('211');plot(1:Nfft,Sxw); xlabel('\omega/\pi');ylabel('e^j^\omega');title('截取信號的頻域');
subplot(212);plot(1:Nfft,abs(Sxw)); xlabel('\omega/\pi');ylabel('|e^j^\omega|');title('截取信號的頻域');
%============================================================================
figure;
subplot(211);
plot(Sxdb); %橫軸1-2000
ylabel('Magnitude (dB)'); xlabel('Frequency (kHz)');title('截取信號的功率譜');
subplot(2,1,2); f1 = (0 : Nfft/2-1)*fs / Nfft / 1000;%取前一半 后一半是翻轉，不必考慮
plot(f1, Sxdb); ylabel('Magnitude (dB)'); xlabel('Frequency (kHz)');title('截取信號的功率譜');
====================================================================================================================================
經典譜估計法1：相關圖法
為了減少譜估計的方差，采用長度為2M-1的窗函數對自相關函數進行截取（聯系上式），得
可使用矩形窗和三角窗。
估計自相關序列：
這里解釋一下，下標之所以是0~L-1 且r關於l=0對稱，是因為數學推導：把-l帶入r(l)依然能得到后面式子的結果。（問的老師）
構成加窗自相關序列：
關於自相關的補充：
計算序列 f(l) 的NFFT（一般選NFFT >2L-1）點DFT/FFT，即為功率譜估計的采樣值:
%進一步處理
r = zeros(2*N/2-1, 1);%(-(N/2-1)~(N/2-1))共2L-1個點 計算自相關%3999*1 %這個不是完全的自相關，只是一半的自相關
for k = 1 : N/2 %從1到2000循環
x1 = x(k : N); %從k到4000
x2 = x(1 : N+1-k); %從1到4001-k
r(N/2+k-1) = x1' * x2 / N;
r(N/2-k+1) = r(N/2+k-1); %r(-k) = r(k)
end
rx = r ;
Sxz1 = fft(rx, Nfft); %DFT
Sxdbz1 = 10*log10(abs(Sxz1(1 : Nfft/2)));%轉換成db
figure;
subplot(211);
plot(1:Nfft,abs(Sxz1)); ylabel('|e^j^\omega|');xlabel('\omega/\pi');title('相關圖法');
subplot(212); plot(f1, Sxdbz1); ylabel('振幅 (dB)'); xlabel('頻率 (kHz)');title('相關圖法（矩形窗)‘）;
%=======================================================================================
w = triang(2*N/2-1)'; %三角窗 加窗后效果
rx = r .* w';
Sxz2 = fft(rx, Nfft);
Sxdbz2 = 10*log10(abs(Sxz2(1 : Nfft/2)));
figure;
plot(f1, Sxdbz2);
ylabel('振幅 (dB)'); xlabel('頻率 (kHz)');
title('相關圖法----加三角窗后');
經典譜估計法2：周期圖法
由本文開始給的定義，功率譜的計算可以是時域信號傅氏變換模平方然后除以時間長度。
但是此方法，當周期圖的方差(當N較大時)，方差：
（可以用多次實驗取平均來緩解） 改進：
- 多個周期圖求平均 把數據記錄切分為K個分段，分別求周期圖，然后求平均。每段長L，偏移量D
（上式@號其實是=號） PA: 周期圖求平均；
Bartlett方法：D=L; Welch方法： D=L/2
Bartlett方法就是把數據分D段，每段fft模平方除以每段長度，再把D段的s相加再平均。
Welch方法就是有重復的分段，具體如下圖：
%----------------周期圖法 Bartlett譜估計--------------%
Sx = zeros(1, Nfft/2);K = 4; L = N/K; %Sx:2000*1 ,L=1000
for k = 1 : K
ks = (k-1)*L + 1; %k=1,ks=1; k=2,ks=1001;
ke = ks + L - 1; %k=1,ke=1000 ;k=2,ke=2000;
X = fft(x(ks:ke), Nfft);
X = (abs(X)).^2; %周期圖法這里要abs + 平方 注意
for i = 1 : Nfft/2 %i=1:2000
Sx(i) = Sx(i) +X(i);
end
end
for i = 1 : Nfft/2
Sx(i) = 10*log10(Sx(i)/(K*L));
end
figure; %
subplot(4,1,1); plot(f1, Sx); ylabel('Magnitude (dB)'); xlabel('Frequency (kHz)'); title('Bartlett Estimate, N=4000, K=4, D=L=1000')
%----------------周期圖法 Welch譜估計--------------%
Nfft = 4000; K = 8; D = fix(Nfft/2 / (K+1));%fix:向0方向靠攏取整 分為K+1格，可以重疊K次做fft %D=222
L = 2*D; %L=444
Sxw = zeros(1, Nfft/2); %1*2000
w = (window('hamming', L))'; %1*444
for k = 1 : K %1*8
ks = (k-1)*D + 1; %k=1,ks=1;k=2,ks=223;k=3,ks=445;k=4,ks=667; k=8,ks=1555
ke = ks + L; %k=1,ke=445;k=2,ke=667 k=8,ke=1999
xk = x(ks:ke)*w; %時域加窗 %k=1,444*1 1*444
X = fft(xk, Nfft);
X = (abs(X)).^2;
for i = 1 : Nfft/2
Sxw(i) = Sxw(i) + X(i); %這里只取前N/2個點 因為后N/2個點是前的翻轉
end
end
for i = 1 : Nfft/2
Sxw(i) = 10*log10(Sxw(i)/(K*L)); %轉換成db
end
figure;
plot(f1, Sxw); ylabel('Magnitude (dB)'); xlabel('Frequency (kHz)'); title('Welch Estimate, N=4000, K=4, D=222, L=444');
語譜圖
語音信號的時頻分布為定義在二維空間的函數，把時頻分布畫成二維灰度圖像的形式，即為語譜圖。
MATLAB 函數 [S, f, t, P] = spectrogram(x, window, noverlap, nfft, fs)； 效果圖
%--------------語譜圖--------------%
bw = 300;
nwin = round(2*fs/bw); %nfft = 512;
nfft = 1024;
xy = filter([1 -0.97], 1, s);
noverlap = nwin - round(length(s) / 500);
% compute and show
figure;
[S, f, t, P] = spectrogram(xy, nwin, noverlap, nfft, fs);
surf(t, f/1000, 10*log10(abs(P)), 'EdgeColor', 'none');
axis xy; axis tight; colormap(jet); view(0, 90); xlabel('Time (s)'); ylabel('Frequency (kHz)');
%title('Broadband Spectrogram');
title('Narrowband Spectrogram');
spectrogram
功能：使用短時傅里葉變換得到信號的頻譜圖。
語法：
[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)
[S,F,T,P]=spectrogram(x,window,noverlap,F,fs)
說明：當使用時無輸出參數，會自動繪制頻譜圖；有輸出參數，則會返回輸入信號的短時傅里葉變
換。當然也可以從函數的返回值S,F,T,P繪制頻譜圖，具體參見例子。
參數：
x---輸入信號的向量。默認情況下，即沒有后續輸入參數，x將被分成8段分別做變換處理，
如果x不能被平分成8段，則會做截斷處理。默認情況下，其他參數的默認值為
window---窗函數，默認為nfft長度的海明窗Hamming
noverlap---每一段的重疊樣本數，默認值是在各段之間產生50%的重疊
nfft---做FFT變換的長度，默認為256和大於每段長度的最小2次冪之間的最大值。
另外，此參數除了使用一個常量外，還可以指定一個頻率向量F
fs---采樣頻率，默認值歸一化頻率
Window---窗函數，如果window為一個整數，x將被分成window段，每段使用Hamming窗函數加窗。
如果window是一個向量，x將被分成length(window)段，每一段使用window向量指定的
窗函數加窗。所以如果想獲取specgram函數的功能，只需指定一個256長度的Hann窗。
Noverlap---各段之間重疊的采樣點數。它必須為一個小於window或length(window)的整數。
其意思為兩個相鄰窗不是尾接着頭的，而是兩個窗有交集，有重疊的部分。
Nfft---計算離散傅里葉變換的點數。它需要為標量。
Fs---采樣頻率Hz，如果指定為[]，默認為1Hz。
S---輸入信號x的短時傅里葉變換。它的每一列包含一個短期局部時間的頻率成分估計，
時間沿列增加，頻率沿行增加。
如果x是長度為Nx的復信號，則S為nfft行k列的復矩陣，其中k取決於window，
如果window為一個標量，則k = fix((Nx-noverlap)/(window-noverlap))
如果window為向量，則k = fix((Nx-noverlap)/(length(window)-noverlap))
對於實信號x，如果nfft為偶數，則S的行數為(nfft/2+1)，如果nfft為奇數，
則行數為(nfft+1)/2，列數同上。
F---在輸入變量中使用F頻率向量，函數會使用Goertzel方法計算在F指定的頻率處計算頻譜圖。
指定的頻率被四舍五入到與信號分辨率相關的最近的DFT容器(bin)中。而在其他的使用nfft
語法中，短時傅里葉變換方法將被使用。對於返回值中的F向量，為四舍五入的頻率，其長度
等於S的行數。
T---頻譜圖計算的時刻點，其長度等於上面定義的k，值為所分各段的中點。
P---能量譜密度PSD(Power Spectral Density)，對於實信號，P是各段PSD的單邊周期估計；
對於復信號，當指定F頻率向量時，P為雙邊PSD。
P矩陣的元素計算公式如下P(I,j)=k|S(I,j)|2，其中的的k是實值標量，定義如下
對於單邊PSD，計算公式如下，其中w(n)表示窗函數，Fs為采樣頻率，在0頻率和奈奎斯特
頻率處，分子上的因子2改為1；
對於雙邊PSD，計算公式如下
如果采樣頻率沒有指定，分母上的Fs由2*pi代替。
spectrogram(...)當調用函數時沒有輸出參數，將會自動繪制各段的PSD估計，繪制的命令如下
surf(T,F,10*log10(abs(P)));
axis tight;
view(0,90);
spectrogram(...,'freqloc')使用freqloc字符串可以控制頻率軸顯示的位置。當freqloc=xaxis
時，頻率軸顯示在x軸上，當freqloc=yaxis時，頻率軸顯示在y軸上，默認是顯示在x軸
上。如果在指定freqloc的同時，又有輸出變量，則freqloc將被忽略。
例.計算並顯示二次掃頻信號的PSD圖，掃頻信號的頻率開始於100Hz，在1s時經過200Hz
T = 0:0.001:2;
X = chirp(T,100,1,200,'q');
spectrogram(X,128,120,128,1E3);
title('Quadratic Chirp');
參考：
https://blog.csdn.net/qq_38559814/article/details/86521602
https://blog.csdn.net/bluehatihati/article/details/84097955
spectrogram函數做短時傅里葉分析：https://blog.csdn.net/zhuguorong11/article/details/77801977