帮我写一段matlab代码基于倒谱法提取语音信号的共振峰

以下是基于倒谱法提取语音信号的共振峰的MATLAB代码：

%% 读取语音信号 [y,fs] = audioread('example.wav');

%% 预处理语音信号 % 消除直流分量 y = y - mean(y);

% 预加重 pre_emphasis = [1 -0.97]; y = filter(pre_emphasis,1,y);

%% 求取倒谱系数 % 分帧 frame_size = 20; frame_shift = 10; nfft = 2^nextpow2(frame_size); num_frames = floor((length(y)-frame_size)/frame_shift) + 1; frames = zeros(nfft,num_frames); for i = 1:num_frames start_idx = (i-1)*frame_shift + 1; end_idx = start_idx + frame_size - 1; frames(:,i) = y(start_idx:end_idx).*hamming(frame_size,nfft); end

% 计算功率谱 power_spectrum = abs(fft(frames,nfft)).^2;

% 计算倒谱系数 cepstrum = ifft(log(power_spectrum)); num_coeffs = 13; % 取13个倒谱系数 cepstrum_coeffs = cepstrum(1:num_coeffs,:);

%% 求取共振峰频率 % LPC分析 lpc_coeffs = zeros(num_coeffs,num_frames); for i = 1:num_frames lpc_coeffs(:,i) = lpc(frames(:,i),num_coeffs-1); end

% 求取极点 [roots_polynomials,~] = roots(lpc_coeffs);

% 保留实部为正的极点 valid_roots = roots_polynomials(imag(roots_polynomials)>=0,:); valid_freqs = acosd(abs(valid_roots));

% 绘制频谱图和共振峰频率 figure; subplot(2,1,1); plot(y); title('语音信号'); xlabel('样本数'); ylabel('幅度');

subplot(2,1,2); freq_scale = (0:nfft/2-1)/nfft*fs; freq_scale = freq_scale(1:num_coeffs); plot(freq_scale,cepstrum_coeffs); hold on; scatter(valid_freqs(:,1),valid_freqs(:,2),'rx'); title('倒谱系数和共振峰频率'); xlabel('频率（Hz）'); ylabel('倒谱系数');