s = 4000; the width of filter windows ts numbers per windows筛选器窗口的宽度 ts每个窗口的数量 wx = ha用中文对这段代码的每一行进行详细分析nningts; generates a hanning filter coefficients产生一个汉宁 滤波器系数 fs = 2000; the sampling rate of
这段代码主要用于计算两个信号之间的相关性和相位差。下面是对每一行代码的详细分析:
-
s = 4000; % the width of filter windows, ts numbers per windows 定义了筛选器窗口的宽度,即每个窗口包含的数据点数。
-
wx = hann(ts); % generates a hanning filter coefficients 生成一个汉宁滤波器系数,用于窗口化数据。
-
fs = 2000; % the sampling rate of probes, [kHz] 探针的采样频率,以kHz为单位。
-
fre_cal = fix((length(data1)-ts)/(ts/2)); % calculate the frequency of analysis 计算信号频率的分析范围。
-
range_butter = [1, 999999]/1000000; % the frequency range of butterworth filter Butterworth滤波器的频率范围。
-
[coef1, coef2] = butter(1, range_butter); % generates the coefficients of butterworth filter 生成Butterworth滤波器的系数。
-
data1_but = filter(coef1, coef2, filter(coef1, coef2, data1)); % data processed by butterworth filter 用Butterworth滤波器处理信号1。
-
data2_but = filter(coef1, coef2, filter(coef1, coef2, data2)); % data processed by butterworth filter 用Butterworth滤波器处理信号2。
-
data1_fft = zeros(fre_cal, ts); % generates a box for the data after FFT 为FFT处理后的数据生成一个盒子。
-
data2_fft = zeros(fre_cal, ts); % generates a box for the data after FFT 为FFT处理后的数据生成一个盒子。
-
disp('系综平均数:') disp(fre_cal) 输出系综平均数。
-
for i = 1 : fre_cal data1_fft(i, :) = data1_but((i-1)*ts/2 + 1 : (i-1)*ts/2 + ts); % read the data for every filter windows data2_fft(i, :) = data2_but((i-1)*ts/2 + 1 : (i-1)*ts/2 + ts); data1_detr(i, :) = detrend(data1_fft(i, :)); % reduce the ripple of data data2_detr(i, :) = detrend(data2_fft(i, :)); data1_fft(i, :) = fft(data1_detr(i, :).*wx', ts); data2_fft(i, :) = fft(data2_detr(i, :).*wx', ts); end 对每个筛选器窗口的数据进行FFT处理,并减去趋势成分。然后将窗口化后的数据进行FFT变换。
-
frequency = (0 : (ts/2 - 1))*(fs/(ts - 2)); % calculate the frequency series. 计算频率序列。
-
mx12 = sum(abs(data1_fft(:, 1:ts/2)).*abs(data2_fft(:, 1:ts/2)))/fre_cal; % calculate (Σ|x1|·|x2|)/fre_cal 计算两个信号之间的相关性。
-
phase = angle(sum(data1_fft(:, 1:ts/2).*conj(data2_fft(:, 1:ts/2))))/pi; % calculate the phase of all filter windows, [rad] 计算所有筛选器窗口的相位差。
-
power_cross = abs(sum(data1_fft(:, 1:ts/2).*conj(data2_fft(:, 1:ts/2))))/fre_cal; % calculate the average cross power of all filter windows 计算所有筛选器窗口的平均交叉功率。
-
power1 = sum(((abs(data1_fft(:, 1:ts/2))).^2))/fre_cal; % calculate the power of data1 计算信号1的功率。
-
power2 = sum(((abs(data2_fft(:, 1:ts/2))).^2))/fre_cal; % calculate the power of data2 计算信号2的功率。
-
corherency = power_cross./mx12; % calculate corherency 计算相关性。
-
freOutput = frequency(2:ts/2); phaseOutput = phase(2 : end); power1Output = power1(2:end); power2Output = power2(2:end); powerCsOutput = power_cross(2:end); CorrOutput = corherency(2:end); 将计算的结果存入输出变量中。
-
figure subplot(511) plot(frequency(2 : ts/2), smooth(power1(2 : end), 3), 'LineWidth', 1), xlim([fBegin, fEnd]) ylabel('power1 (a.u)') subplot(512) plot(frequency(2 : ts/2), smooth(power2(2 : end), 3), 'LineWidth', 1), xlim([fBegin, fEnd]) ylabel('power2 (a.u)') subplot(513) plot(frequency(2 : ts/2), power_cross(2 : end), 'LineWidth', 1), xlim([fBegin, fEnd]) ylabel('cross power (a.u)') subplot(514) plot(frequency(2 : ts/2), smooth(corherency(2 : end), 5), 'LineWidth', 1), xlim([fBegin, fEnd]) ylabel('corherency') subplot(515) plot(frequency(2 : ts/2), smooth(phase(2 : end), 3), 'LineWidth', 1), xlim([fBegin, fEnd]) hold on plot([0, 500], [0, 0], 'b--', 'linewidth', 1), ylim([-1, 1]) ylabel('phase (rad/\pi)') xlabel('f (kHz)') 绘制相关性和相位差的图像。subplot函数将不同的图像绘制在一个图中。每个图像的横轴为频率,纵轴为功率、相关性或相位差。
原文地址: https://www.cveoy.top/t/topic/bpOj 著作权归作者所有。请勿转载和采集!