MATLAB音频信号抽样与恢复代码解析

% 清除命令窗口和工作区变量 clc clear all % 读取音频信号 [x, fs] = audioread('D:\USERS\Desktop\音频处理\om.mp3'); % 确定限带频率 % 计算FFT和功率谱密度（PSD） N = length(x); X = fft(x)/N; PSD = abs(X).^2/N; % 确定限带频率 freq = linspace(0,fs/2,floor(N/2)+1); % 创建一个频率向量 P = PSD(1:length(freq)); % 提取PSD中的有限频率 [~,ind] = sort(P,'descend'); % 按降序排列功率谱密度 cumP = cumsum(P(ind)); % 计算排序后功率谱密度的累积分布 cumP = cumP/cumP(end); % 归一化 threshold = 0.9; % 限定阈值为0.9 f_low = freq(ind(find(cumP>=threshold,1))); % 获取低频率限制 f_high = freq(ind(find(cumP>=threshold,1,'last'))); % 获取高频率限制 disp(['Signal bandwidth is between ', num2str(f_low), ' Hz and ', num2str(f_high), ' Hz.']); % 输出限带频率信息 % 确定信号最高频率 f = (0:N-1)(fs/N); % f是每个采样点的频率 [maxValue, maxIndex] = max(abs(X)); % 获取FFT中的最大幅值和相应的下标 freq = f(maxIndex); % 根据下表获取信号最高频率 fc=2freq(1,1); % 抽样频率需大于信号最高频率的两倍 fc_r=round(fc); % 取整数 numStr = num2str(fc_r); % 取出第一位数字 firstDigit = str2double(numStr(1)); % 如果第一位数字是1-9之间的数，就将整数向上取整到这个数字的最高位 if firstDigit > 0 && firstDigit <= 9 rounded_fc = ceil(fc_r / 10^(length(numStr)-1)) * 10^(length(numStr)-1); else % 如果第一位数字是0，则向上取整到第一个非零数字的最高位 for i = 2:length(numStr) if str2double(numStr(i)) > 0 rounded_fc = ceil(fc_r/ 10^(length(numStr)-i)) * 10^(length(numStr)-i); return; end end rounded_fc = 0; % 如果整数是0，则无法取整 end % 三种抽样频率 Fs_1 = rounded_fc; % 等于两倍奈奎斯特抽样频率 Fs_2 = 2rounded_fc; % 大于两倍奈奎斯特抽样频率 Fs_3 = 0.8rounded_fc; % 小于两倍奈奎斯特抽样频率 % 抽样 n_1 = length(x); % 两倍抽样频率抽样的点数 t_1 = (0:n_1-1)/fs; n_2 = round(n_1 * Fs_2 / Fs_1); % 大于两倍抽样频率抽样的点数 t_2 = (0:n_2-1)/Fs_2; % 抽样信号的时间向量 n_3 = round(n_1 * Fs_3 / Fs_1); % 小于两倍抽样频率抽样的点数 t_3 = (0:n_3-1)/Fs_3; % 抽样信号的时间向量 x_1 = x; % 初始化抽样信号1 x_2 = zeros(n_2, 1); % 初始化抽样信号2 x_3 = zeros(n_3, 1); % 初始化抽样信号3 % 生成以Fs_2抽样的信号 for i = 1:n_2 x_2(i) = x(ceil(i * Fs_1 / Fs_2)); end % 生成以Fs_3抽样的信号 for i = 1:n_3 x_3(i) = x(ceil(i * Fs_1 / Fs_3)); end % 恢复 y_1 = x_1; % 原始信号 y_2 = zeros(n_1, 1); % 初始化恢复信号2 y_3 = zeros(n_1, 1); % 初始化恢复信号3 % 以Fs_2抽样的信号恢复 for i = 1:n_2 y_2(ceil(i * Fs_1 / Fs_2)) = x_2(i); end % 以Fs_3抽样的信号恢复 for i = 1:n_3 y_3(ceil(i * Fs_1 / Fs_3)) = x_3(i); end % 时域图 figure(1) subplot(3,2,1); plot(t_1, x_1); title(sprintf('Sampled Signal (Fs = %dHz)',Fs_1)); xlabel('Time (s)'); ylabel('Amplitude'); subplot(3,2,2); plot(t_1, y_1); title(sprintf('Recovered Signal (Fs = %dHz)',Fs_1)); xlabel('Time (s)'); ylabel('Amplitude'); subplot(3,2,3); plot(t_2, x_2); title(sprintf('Sampled Signal (Fs = %dHz)', Fs_2)); xlabel('Time (s)'); ylabel('Amplitude'); subplot(3,2,4); plot(t_1, y_2); title(sprintf('Recovered Signal (Fs = %dHz)', Fs_2)); xlabel('Time (s)'); ylabel('Amplitude'); subplot(3,2,5); plot(t_3, x_3); title(sprintf('Sampled Signal (Fs = %dHz)', Fs_3)); xlabel('Time (s)'); ylabel('Amplitude'); subplot(3,2,6); plot(t_1, y_3); title(sprintf('Recovered Signal (Fs = %dHz)', Fs_3)); xlabel('Time (s)'); ylabel('Amplitude'); % 频谱图 N = length(x_1); f = fs*(0:(N/2))/N; X_1 = fft(x_1)/N; X_2 = fft(x_2)/N; X_3 = fft(x_3)/N; Y_1 = fft(y_1)/N; Y_2 = fft(y_2)/N; Y_3 = fft(y_3)/N; figure(2) subplot(3,2,1); plot(f, 2abs(X_1(1:N/2+1))); title(sprintf('Sampled Signal (Fs = %dHz)', Fs_1)); xlabel('Frequency (Hz)'); ylabel('Magnitude'); subplot(3,2,2); plot(f, 2abs(Y_1(1:N/2+1))); title(sprintf('Recovered Signal (Fs = %dHz)', Fs_1)); xlabel('Frequency (Hz)'); ylabel('Magnitude'); subplot(3,2,3); plot(f, 2abs(X_2(1:N/2+1))); title(sprintf('Sampled Signal (Fs = %dHz)', Fs_2)); xlabel('Frequency (Hz)'); ylabel('Magnitude'); subplot(3,2,4); plot(f, 2abs(Y_2(1:N/2+1))); title(sprintf('Recovered Signal (Fs = %dHz)', Fs_2)); xlabel('Frequency (Hz)'); ylabel('Magnitude'); subplot(3,2,5); plot(f, 2abs(X_3(1:N/2+1))); title(sprintf('Sampled Signal (Fs = %dHz)', Fs_3)); xlabel('Frequency (Hz)'); ylabel('Magnitude'); subplot(3,2,6); plot(f, 2abs(Y_3(1:N/2+1))); title(sprintf('Recovered Signal (Fs = %dHz)', Fs_3)); xlabel('Frequency (Hz)'); ylabel('Magnitude'); figure(3) % 原始信号 plot(x, x_1); title(sprintf('Sampled Signal (Fs = %dHz)',Fs_1)); xlabel('Time (s)'); ylabel('Amplitude');