使用matlab对信号进行经典谱估计

部分内容摘自 https://blog.csdn.net/liusandian/article/details/53386581

功率谱和频谱

先简要说计算：
功率谱：信号先自相关再作FFT
频 谱：信号直接作FFT

功率谱：
信号的传播都是看不见的，但是它以波的形式存在着，这类信号会产生功率，单位频带的信号功率就被称之为功率谱。它可以显示在一定的区域中信号功率随着频率变化的分布情况。

功率谱可以从两方面来定义：
一个 是自相关函数的傅立叶变换；（维纳辛钦定理）
另一个 是时域信号傅氏变换模平方然后除以时间长度。（来自能量谱密度）
根据parseval定理，信号傅氏变换模平方被定义为能量谱，能量谱密度在时间上平均就得到了功率谱。

频谱：
频谱是常常指信号的Fourier变换。

（1-7 作者：Yorkxu）转载的理解：
（1）信号通常分为两类：能量信号和功率信号；

• 能量信号：又称能量有限信号，是指在所有时间上总能量不为零且有限的信号。
• 功率信号：它的能量为无限大，它对通信系统的性能有很大影响，决定了无线系统中发射机的电压和电磁场强度。

（2）一般来讲，能量信号其傅氏变换收敛（即存在），而功率信号傅氏变换通常不收敛（当然，若信号存在周期性，可引入特殊数学函数（Delta）表征傅氏变换的这种非收敛性）；

（3）信号是信息的搭载工具，而信息与随机性紧密相关，所以实际信号多为随机信号，这类信号的特点是状态随机性随时间无限延伸，其样本能量无限。换句话说，随机信号（样本）大多属于功率信号而非能量信号，它并不存在傅氏变换，亦即不存在频谱；

（4）若撇开搭载信息的有用与否，随机信号又称随机过程，很多噪声属于特殊的随机过程，它们的某些统计特性具有平稳性，其均值和自相关函数具有平稳性。对于这样的随机过程，自相关函数蜕化为一维确定函数，前人证明该确定相关函数存在傅氏变换；

（5）能量信号频谱通常既含有幅度也含有相位信息；幅度谱的平方（二次量纲）又叫能量谱（密度），它描述了信号能量的频域分布；功率信号的功率谱（密度）描述了信号功率随频率的分布特点（密度：单位频率上的功率），业已证明，平稳信号功率谱密度恰好是其自相关函数的傅氏变换。对于非平稳信号，其自相关函数的时间平均（对时间积分，随时变性消失而再次退变成一维函数）与功率谱密度仍是傅氏变换对；

（6）实际中我们获得的往往仅仅是信号的一段支撑，此时即使信号为功率信号，截断之后其傅氏变换收敛，但此变换结果严格来讲不属于任何“谱”（进一步分析可知它是样本真实频谱的平滑：卷积谱）；

（7）对于（6）中所述变换若取其幅度平方，可作为平稳信号功率谱（密度）的近似，是为经典的“周期图法”；

补充：（8） 一个信号的频谱,只是这个信号从时域表示转变为频域表示,只是同一种信号的不同的表示方式而已；而功率谱是从能量的观点对信号进行的研究,其实频谱和功率谱的关系归根揭底还是信号和功率,能量等之间的关系。

谱估计

功率谱估计一般分成两大类：
经典谱估计，也称为非参数谱估计。
现代谱估计，也称为参数谱估计。

根据维纳-辛钦定理，平稳随机过程的功率谱等于自相关序列的DTFT:

matlab做信号预处理

一开始博主录音了内容为“现代信号处理”的话音，并存为单声道wav做处理，但是因为信号过长 且发音不标准 导致实验效果受影响，所以减少话音长度，仅仅读取“代”的话音作为分析样本。
上源码：

close all;clear all;
% read speech waveform from a file
[s, fs] = audioread('代.wav');

% set analysis parameters, pre-emphasise and windowing
%根据话音位置，取8192可以把话音主要部分加入其中
N = 8192;
Nfft = 8192;
n0 = 1000;
x = s(n0 : n0+N-1);
x1 = filter([1 -0.97], 1, x); %预加重 滤波器
%作用：消除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(4,1,1);
plot(s);  xlim([0 length(s)]); ylim([-0.65 0.65]);
ylabel('Amplitude');  xlabel('Time (n)');
subplot(4,1,2);
plot(xw); xlim([0 length(x)]); ylim([-0.225 0.225]);
ylabel('Amplitude');  xlabel('Time (n)'); 

f = (0 : Nfft/2-1)*fs / Nfft / 1000;%取前一半 后一半是翻转
subplot(2,1,2);
plot(f, Sxdb);
ylabel('Magnitude (dB)');  xlabel('Frequency (kHz)');

效果图：

从上至下：原始话音、经过预处理后的话音、话音的理论上的功率谱

经典谱估计法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个点 计算自相关
for k = 1 : N/2
    x1 = x(k : N);
    x2 = x(1 : N+1-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);
Sxdbz1 = 10*log10(abs(Sxz1(1 : Nfft/2)));%转换成db

subplot(4,1,3);
plot(f, Sxdbz1);
ylabel('Magnitude (dB)');  xlabel('Frequency (kHz)');
title('Correlogram (Rectangle)');

%----------------相关图法 三角窗 --------------%
w = triang(2*N/2-1)'; %三角窗 加窗后效果
rx = r .* w';
Sxz2 = fft(rx, Nfft);
Sxdbz2 = 10*log10(abs(Sxz2(1 : Nfft/2)));
subplot(4,1,4);
plot(f, Sxdbz2);
ylabel('Magnitude (dB)');  xlabel('Frequency (kHz)');
 title('Correlogram (Triangular)');

效果图：
下两幅是相关图法的结果，分别使用了矩形窗和三角窗

经典谱估计法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;
for k = 1 : K
    ks = (k-1)*L + 1;
    ke = ks + L - 1;
    X = fft(x(ks:ke), Nfft);
    X = (abs(X)).^2;          %周期图法这里要abs + 平方 注意
    for i = 1 : Nfft/2
        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(f, Sx);
ylabel('Magnitude (dB)');  xlabel('Frequency (kHz)');
title('Bartlett Estimate, N=4096, K=4, D=L=1024')

%----------------周期图法 Welch谱估计--------------%
Nfft = 8192; 
K = 4; 
D = fix(Nfft/2 / (K+1));%向0方向靠拢取整 分为K+1格，可以重叠K次做fft
L = 2*D;
Sxw = zeros(1, Nfft/2);
w = (window('hamming', L))';
for k = 1 : K
    ks = (k-1)*D + 1;
    ke = ks + L - 1;
    xk = x(ks:ke) .* w; %时域加窗
    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
subplot(4,1,2);  plot(f, Sxw);
ylabel('Magnitude (dB)'); xlabel('Frequency (kHz)');
title('Welch Estimate, N=4096, K=4, D=819, L=1638');

效果图是效果图的第一和第二张。

语谱图

语音信号的时频分布为定义在二维空间的函数，把时频分布画成二维灰度图像的形式，即为语谱图。

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');

END

致敬北邮尹霄丽老师！公式和很多概念都引用老师现代信号处理课程的教案，老师上课讲的很细致很明白，很荣幸能作为老师的学生。