当前位置：首页 > news >正文

PNCC（Power-Normalized Cepstral Coefficients）— MATLAB 实现

news 2026/6/16 16:08:04

对应论文：C. Kim & R. M. Stern, Power-Normalized Cepstral Coefficients (PNCC) for Robust Speech Recognition, IEEE/ACM TASLP 2010.

一、PNCC 流水线到底做了什么

语音x[n]│├─ ① 预加重  y[n] = x[n] - 0.97·x[n-1]│├─ ② 分帧  25.6ms Hamming / 10ms hop / NFFT=1024 @16kHz│         → 短时功率谱 |X[m,k]|²│├─ ③ Mel/Gammatone 滤波器组（论文用40通道）→ 各通道功率 P[m,l]│                                         l=1..L (L=40)│├─ ④ ★Medium-Time Power（5帧滑动均值）│       Q̃[m,l] = (1/5)Σ_{i=-2}^{+2} P[m+i, l]│├─ ⑤ ★Asymmetric Noise Suppression (ANS)│       → 不对称低通滤出噪声底 Q̃_le│       → R[m,l] = max( Q̃[m,l] - c·Q̃_le , T[m,l] )│       → Temporal Masking（峰值跟踪）│├─ ⑥ ★Time-Frequency / Mean Power Normalization│       → 除以帧总功率均值（或平滑transfer fn）│├─ ⑦ Power-Law 非线性  (·)^(1/15)     ← 替换了 MFCC 的 log(·)│├─ ⑧ DCT → c0..c12└─ ⑨ CepLifter → 输出 13 维 PNCC（可追加 Δ/ΔΔ）

二、主函数 `extract_pncc.m`

function [pncc, aux] = extract_pncc(x, fs, varargin)
% -------------------------------------------------------------
% PNCC特征提取（手写，无工具箱）
%   [pncc, aux] = extract_pncc(x, fs)
%   [pncc, aux] = extract_pncc(x, fs, 'numceps',13, 'do_delta',true)
%
% 输入:
%   x  : 语音信号（double列向量，任意幅值）
%   fs : 采样率（推荐16000）
% 输出:
%   pncc : [Nframes × numceps]  静态系数（通常13维）
%   aux  : 中间量结构体（可可视化调试）
% -------------------------------------------------------------%% ---------- 参数 ----------
p.numceps   = 13;      % 输出cepstral维数（含c0=第1列）
p.preemph   = 0.97;    % 预加重系数
p.win_ms    = 25.6;    % 帧长 ms（论文25.6ms @16kHz）
p.hop_ms    = 10;      % 帧移 ms
p.Nfft      = 1024;    % FFT点数
p.L         = 40;      % 滤波器组通道数（40）
p.f_low     = 200;     % 滤波器组最低频率 Hz
p.f_high    = 8000;    % 最高频率 Hz（@16k常设为fs/2）
p.power_exp = 1/15;    % 幂律指数论文用 1/15
p.lifter    = 22;       % cep lifter系数（可选）
p.do_delta  = false;    % 是否追加Δ/ΔΔ → 39维
p.noise_c   = 0.5;     % ANS公式里 c（论文建议0~1之间，常取0.5）
p.floor_eps = 1e-12;    % 数值地板% 解析varargin
for i=1:2:numel(varargin)if isfield(p, varargin{i}), p.(varargin{i}) = varargin{i+1}; end
end%% ---------- ① 预加重 ----------
x = x(:);
x = x - p.preemph * [0; x(1:end-1)];   % y[n]-0.97y[n-1]
x = x * 0.99 / max(abs(x), 1e-12);    % 粗归一化（不影响特征比）%% ---------- ② 分帧 ----------
wlen = round(p.win_ms/1000 * fs);
hop  = round(p.hop_ms/1000 * fs);
Nfft = p.Nfft;
nwin = floor((length(x)-wlen)/hop) + 1;% Hamming窗
win = 0.54 - 0.46*cos(2*pi*(0:wlen-1)'/(wlen-1));frames = zeros(nwin, wlen);
for m = 1:nwinframes(m,:) = x((m-1)*hop + (1:wlen)) .* win.';
end%% ---------- ③ 短时功率谱 ----------
F = Nfft/2 + 1;
Pxx = zeros(nwin, F);      % |X|²
for m = 1:nwinX = fft(frames(m,:), Nfft);X = X(1:F);Pxx(m,:) = real(X).^2 + imag(X).^2;
end%% ---------- ③b Mel三角滤波器组 ----------
fbank = build_triangular_mel_fbank(fs, Nfft, p.L, p.f_low, p.f_high);
% fbank : [F × L]  , 每列是一个mel通道的权向量（已归一化可选）
% 通道功率: P[m,l] = Σ_k Pxx[m,k] * fbank[k,l]
P_ch = zeros(nwin, p.L);
for m = 1:nwinP_ch(m,:) = (Pxx(m,:) * fbank).';   % 1×L
end%% ---------- ④ Medium-Time Power ----------
M = 2;   % 5帧窗口：m-2..m+2
Q_tilde = zeros(nwin, p.L);
for m = 1:nwinidx = max(1, m-M) : min(nwin, m+M);Q_tilde(m,:) = mean(P_ch(idx,:), 1);
end%% ---------- ⑤ Asymmetric Noise Suppression (ANS) ----------
% Step-1: 不对称低通（一阶IIR：α≈0.999, β=0.5 论文写法）
% 这里用因果的一阶IIR近似 Q_le[m] = α·Q_le[m-1] + (1-α)·Q_tilde[m]
alpha_asym = 0.999;
beta_asym  = 0.5;
Q_tilde_le = zeros(nwin, p.L);
Q_tilde_le(1,:) = Q_tilde(1,:);
for m = 2:nwinQ_tilde_le(m,:) = alpha_asym * Q_tilde_le(m-1,:) + (1-alpha_asym) * Q_tilde(m,:);
end
% 可选：再做一次反向（非因果）使更接近论文"asymmetric lowpass"
Q_le_rev = zeros(nwin, p.L);
Q_le_rev(end,:) = Q_tilde_le(end,:);
for m = nwin-1:-1:1Q_le_rev(m,:) = alpha_asym * Q_le_rev(m+1,:) + (1-alpha_asym) * Q_tilde_le(m,:);
end
Q_le = (Q_tilde_le + Q_le_rev)/2;% Step-2:  subtract + half-wave rectifier
% R[m,l] = max( Q_tilde[m,l] - c·Q_le[m,l] ,  ε )  （论文T是小的floor）
c = p.noise_c;
R = Q_tilde - c .* Q_le;
R(R < p.floor_eps) = p.floor_eps;   % half-wave rectifier + floor% Step-3: Temporal Masking（论文：每通道找局部峰值，低于峰值衰减）
R_masked = zeros(size(R));
half_win = 1;  % 3帧掩蔽窗口（可调）
for l = 1:p.Lfor m = 1:nwin% 局部峰值（含自己）idx = max(1,m-half_win):min(nwin,m+half_win);peak = max(R(idx,l));if R(m,l) >= 0.8*peakR_masked(m,l) = R(m,l);elseR_masked(m,l) = p.floor_eps;   % 被mask掉的区域只留floorendend
end%% ---------- ⑥ Mean Power Normalization ----------
% 论文：每帧总功率 S[m] = Σ_l R_masked[m,l]，再除 S[m] 或 S^a
S = sum(R_masked, 2); S(S<1e-12)=1e-12;
R_norm = R_masked ./ S;%% ---------- ⑦ Power-Law Nonlinearity ----------
% 论文用 exponent = 1/15 ≈ 0.0667（不是log！）
R_pow = (R_norm + p.floor_eps).^p.power_exp;% 到这一步已经很像"spectrogram"了，可以可视化
melE_final = R_pow;%% ---------- ⑧ DCT → Cepstra ----------
% 手写DCT-II（正交）
cep = zeros(nwin, p.L);
for m = 1:nwinvec = melE_final(m,:).';cep(m,:) = dct_ii(vec(:), p.L).';   % 1×L
endpncc = cep(:,1:min(p.numceps,p.L));%% ---------- ⑨ CepLifter ----------
if p.lifter > 0lifter_vec = 1 + 0.5*p.lifter*sin(pi*(0:size(pncc,2)-1)/p.lifter);pncc = pncc .* lifter_vec;
end%% ---------- Δ / ΔΔ ----------
if p.do_deltaDelta = diff(pncc,1,1); Delta = [Delta(1,:); Delta];Delta2 = diff(Delta,1,1); Delta2 = [Delta2(1,:); Delta2];pncc = [pncc, Delta, Delta2];  % → 39维
end%% ---------- 打包aux（可视化/调试）----------
aux.Pxx = Pxx; aux.fbank = fbank; aux.P_ch = P_ch;
aux.Q_tilde = Q_tilde; aux.Q_le = Q_le; aux.R_after_ANS = R;
aux.R_masked = R_masked; aux.R_norm = R_norm; aux.melE = melE_final;
aux.frames = frames; aux.wlen = wlen; aux.hop = hop;
end

三、必须的辅助函数

3.1 三角 Mel 滤波器组

function fbank = build_triangular_mel_fbank(fs, Nfft, L, f_low, f_high)
% fbank: [(Nfft/2+1) × L]
% 用 mel 刻度均匀分布中心频率
mel = @(f) 2595*log10(1 + f/700);
mel_low = mel(f_low); mel_high = mel(f_high);
mel_c = linspace(mel_low, mel_high, L+2);   % L+2个端点
hz_c = 700*(10.^(mel_c/2595)-1);           % 转回Hzbin_cntr = (0:(Nfft/2)) * fs/Nfft;           % FFT bin中心频率(Hz)fbank = zeros(length(bin_cntr), L);
for l = 1:Lleft  = hz_c(l);cent  = hz_c(l+1);right = hz_c(l+2);for k = 1:length(bin_cntr)f = bin_cntr(k);if f>=left && f<=centfbank(k,l) = (f-left)/(cent-left);elseif f>cent && f<=rightfbank(k,l) = (right-f)/(right-cent);endend% 可选：面积归一化（让每个通道能量守恒）% fbank(:,l) = fbank(:,l) / max(sum(fbank(:,l)), 1e-12);
end
end

3.2 手写 DCT-II（正交）

function C = dct_ii(x, N)
% x: N×1
C = zeros(N,1);
for k = 0:N-1ck = (k==0)*1/sqrt(N) + (k>0)*sqrt(2/N);idx = (0:N-1).';C(k+1) = ck * sum( x .* cos( pi*(2*idx+1)*k/(2*N) ));
end
end

四、Demo

%% demo_pncc.m
clear; clc; close all;% ---- 1) 造一段带噪语音（或换 audioread('your.wav')）
fs = 16000;
t = 0:1/fs:2;
sp = 0.3*sin(2*pi*200*t) .* (1+0.5*sin(2*pi*3*t)) + ...0.15*sin(2*pi*800*t);
noise = 0.06*randn(size(t));
x = sp + noise;
x = x(:);% ---- 2) 提取PNCC ----
[pncc, aux] = extract_pncc(x, fs, 'numceps',13, 'do_delta',false);fprintf('PNCC size = %d × %d\n', size(pncc,1), size(pncc,2));
fprintf('帧0  cep = [%s]\n', sprintf(' %.4f', pncc(1,:)));% ---- 3) 画关键图 ----
figure('Color','w','Position',[100 80 1200 700]);subplot(3,3,1); imagesc(log10(aux.Pxx+1e-12)); axis xy;xlabel('FFT bin'); ylabel('frame'); title('log|Pxx|'); colorbar;subplot(3,3,2); imagesc(log10(aux.P_ch+1e-12)); axis xy;xlabel('channel l'); ylabel('frame'); title('通道功率 P[m,l]'); colorbar;subplot(3,3,3); imagesc(log10(aux.R_after_ANS+1e-12)); axis xy;xlabel('channel l'); ylabel('frame'); title('ANS后 R[m,l]'); colorbar;subplot(3,3,4); imagesc(aux.R_norm); axis xy;xlabel('channel l'); ylabel('frame'); title('Mean-Power-Norm R_norm'); colorbar;subplot(3,3,5); plot(pncc(:,1),'b','LineWidth',1.5);xlabel('frame'); ylabel('c0'); title('c0（0th cepstral）'); grid on;subplot(3,3,6); imagesc(pncc(:,1:13)); axis xy;xlabel('cep coeff'); ylabel('frame'); title('PNCC 特征矩阵'); colorbar;sgtitle('PNCC 全流程（Kim & Stern 2010 手写 MATLAB）');

参考代码 PNCC特征值提取 www.youwenfan.com/contentcnv/81293.html

五、和 MFCC 的核心区别

模块	MFCC	PNCC
滤波器组	三角 Mel	三角 Mel（论文推荐 Gammatone）
非线性	log(·)	(·)^(1/15)（幂律，贴合听神经 rate-intensity）
多帧处理	无（纯短时）	Medium-time（5帧）→ ANS噪声抑制 → 时频归一化
对平稳噪声	靠log压制	靠不对称估计噪声底做减法，更不容易把语音当噪声消掉