功率譜和頻譜：函數
功率譜：信號自相關後FFT工具
頻譜：信號直接FFT性能
功率譜：spa
信號的傳播都是看不見的，可是它以波的形式存在着，這類信號會產生功率，單位頻帶的信號功率就被稱之爲功率譜。它能夠顯示在必定的區域中信號功率隨着頻率變化的分佈狀況。.net
功率譜能夠從兩方面來定義：3d
一個 是自相關函數的傅立葉變換；（維納辛欽定理）orm
另外一個 是時域信號傅氏變換模平方而後除以時間長度。（來自能量譜密度） 根據parseval定理，信號傅氏變換模平方被定義爲能量譜，能量譜密度在時間上平均就獲得了功率譜。 頻譜：blog
頻譜是經常指信號的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