用matlab编写一段计算THD的程序

在MATLAB中，全称THD（Total Harmonic Distortion）是指非线性信号的总谐波失真度，通常用于评估交流电源或信号质量。以下是一个简单的MATLAB函数，用于计算单通道信号的THD：

function thd = calculate_THD(sinusoidal_signal, distorted_signal)
% sinusoidal_signal: 清洁正弦波参考信号
% distorted_signal: 受到干扰的信号

% 首先，我们需要对两个信号进行归一化处理
clean_norm = max(abs(sinusoidal_signal));
distorted_norm = max(abs(distorted_signal));

% 计算归一化的信号
sinusoidal_normalized = sinusoidal_signal ./ clean_norm;
distorted_normalized = distorted_signal ./ distorted_norm;

% 求解FFT（快速傅立叶变换）并取绝对值
fft_sinusoidal = abs(fft(sinusoidal_normalized));
fft_distorted = abs(fft(distorted_normalized));

% 确定频率轴
fs = 1000; % 假设采样率
frequencies = (0:length(fft_sinusoidal) - 1) * fs / length(fft_sinusoidal);

% 截取只保留第一半（奇次谐波）的频谱，因为偶数次谐波是零分量
half_length = floor(length(fft_sinusoidal)/2);
harmonics = frequencies(1:half_length);

% 计算每个谐波的幅度，并找到最大奇数次谐波
harmonic_amplitudes = fft_sinusoidal(1:half_length);
harmonic_index = find(harmonic_amplitudes == max(harmonic_amplitudes));

% 计算总谐波失真
total_harmonic_power = sum(harmonic_amplitudes(harmonic_index));
fundamental_power = fft_sinusoidal(end); % 基准周期功率
thd = total_harmonic_power / fundamental_power;

% 返回THD
end

使用这个函数的例子：

% 假设你有清洁正弦波和受干扰信号的数据
clean_wave = sin(2*pi*50*[0:0.001:1]);
distorted_wave = clean_wave + 0.1*sin(2*pi*100*[0:0.001:1]);

% 调用函数计算THD
thd_value = calculate_THD(clean_wave, distorted_wave);
disp(['The THD of the signal is: ', num2str(thd_value)]);