freq_responce/scripts/getResponse.js

218 lines
11 KiB
JavaScript
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

export function getResponse(args) {
let force = args.spec_force; //спектр силы
let vibro = args.spec_vibro; //спектр вибрации
let type = args.type; //способ вычисления отклика
let view = args.view; //способ отображения спектра в линейных единицах или дБ
/* Алгоритм
1. Получить спектры ускорения и силы в комплексном виде (вибропреобразователь и молоток).
2. Преобразовать спектр силы (молотка) в комплексно-сопряжённый вид.
3. Получить кросс-спектр (перекрёстный) между ускорением и силой.
4. Вычислить деление кросс-спектра на спектр силы.
*/
/* Кросс-спектры
Сигнал силы в комплексно-сопряжённом виде
сигнал вибрации умножается на сигнал силы, а НЕ наоборот
умножение комплексных чисел по формуле:
(a1+b1i)*(a2+b2i)=(a1*a2-b1*b2)*(a1*b2+b1*a2)i
(a1*a2-b1*b2) - действительная часть
(a1*b2+b1*a2)i - мнимая часть
*/
let n = force.data.length; // длина массива спектра
let cross_real = new Float64Array(n); // типизированный массив для действительной части
let cross_imag = new Float64Array(n); // типизированный массив для мнимой части
for (let i = 0; i < n; i++) {
const vReal = vibro.real[i]; // действительная часть спектра вибрации
const vImag = vibro.imag[i]; // мнимая часть спектра вибрации
const fReal = force.real[i]; // действительная часть спектра силы
const fImag = -force.imag[i]; // мнимая часть спектра силы сразу берём с обратным знаком
cross_real[i] = vReal * fReal - vImag * fImag;
cross_imag[i] = vReal * fImag + vImag * fReal;
}
// let cross_real2 = getCrossSpectrum(vibro.real, vibro.imag, force.real, force.imag).reXY;
// gtl.log.info('Амплитуда кросс', cross_real)
// gtl.log.info('Амплитуда кросс 2', cross_real2)
/* Вычисление частотного отклика
деление комплексных чисел по формуле:
(a+bi)/(c+di)=(a*c-b*d)*(a*d+b*c)i
(a*c-b*d) - действительная часть
(a*d+b*c)i - мнимая часть
*/
let resp_real = new Float64Array(n); // типизированный массив для действительной части
let resp_imag = new Float64Array(n); // типизированный массив для мнимой части
if (type == 0) {
// Вариант 1: Частотный отклик как деление кросс-спектра вибрации и силы на автоспектр силы
// ОШИБКА в моей формуле вычисления частотного отклика
// ИПРАВИЛ на:
// деление комплексных чисел по формуле:
// (a+bi)/(c+di)=(a*c + b*d)/(c^2 + d^2) + (c*b - a*d)/(c^2 + d^2)i
// (a*c + b*d)/(c^2 + d^2) - действительная часть
// (c*b - a*d)/(c^2 + d^2)i - мнимая часть
for (let i = 0; i < n; i++) {
const cReal = cross_real[i]; // действительная часть кросс спектра
const cImag = cross_imag[i]; // мнимая часть кросс спектра
const fReal = force.real[i]; // действительная часть спектра силы
const fImag = force.imag[i]; // мнимая часть спектра силы
const denominator = fReal * fReal + fImag * fImag;
resp_real[i] = (cReal * fReal + cImag * fImag) / denominator;
resp_imag[i] = (fReal * cImag - cReal * fImag) / denominator;
}
} else {
// Вариант 2: Частотный отклик как деление автоспектра вибрации на автоспектр силы
// Канал[j] / Канал[0]
for (let i = 0; i < n; i++) {
const vReal = vibro.real[i]; // действительная часть спектра вибрации
const vImag = vibro.imag[i]; // мнимая часть спектра вибрации
const fReal = force.real[i]; // действительная часть спектра силы
const fImag = force.imag[i]; // мнимая часть спектра силы
const denominator = fReal * fReal + fImag * fImag;
resp_real[i] = (vReal * fReal + vImag * fImag) / denominator;
resp_imag[i] = (fReal * vImag - vReal * fImag) / denominator;
}
}
// запись массива амплитуды и фазы
let resp_ampl = new Float64Array(n); // типизированный массив для амплитуды
let resp_phase = new Float64Array(n); // типизированный массив для фазы
const RAD_TO_DEG = 180 / Math.PI; //перевод радиан в градусы
for (let i = 0; i < n; i++) {
const rReal = resp_real[i]; // действительная часть спектра отклика
const rImag = resp_imag[i]; // мнимая часто спектра отклика
resp_ampl[i] = Math.sqrt(rReal * rReal + rImag * rImag);
if (view > 0) { resp_ampl = resp_ampl.map((item) => (item = 20 * Math.log10(item / 1e-6))) }; //переводим в дБ
resp_phase[i] = Math.atan2(rImag, rReal) * RAD_TO_DEG;
}
// Приводим фазу к диапазону [0°, 360°]
for (let i = 0; i < n; i++) {
resp_phase[i] = resp_phase[i] % 360;
if (resp_phase[i] < 0) { resp_phase[i] += 360 }
}
let resp2 = getFrequencyResponse(vibro.real, vibro.imag, force.real, force.imag);
let resp_ampl2 = resp2.magnitude;
let resp_phase2 = resp2.phase;
// gtl.log.info('Амплитуда отклика', resp_ampl[0]);
// gtl.log.info('Амплитуда отклика 2', resp_ampl2[0]);
// gtl.log.info('Фаза отклика', resp_phase[0]);
// gtl.log.info('Фаза отклика 2', resp_phase2[0]);
return {
data: resp_ampl2,
phase: resp_phase2,
resolution: force.resolution
}
}; //рассчет перекрестного спектра
export function getAvgArray(arrays, remove) {
if (!arrays || arrays.length === 0) return []; //защита от пустого входного массива
const length = arrays[0].length; //длина вложенного массива
const result = new Array(length); //новый массив фиксированной длины, чтобы не пушить
const newArrays = arrays.filter((_, i) => !remove.includes(i)); //исключаем массивы по указанным индексам
const numArrays = newArrays.length; //количество вложенных массивов
if (newArrays.length === 0) return []; //защита если после фильтрации ничего не осталось
for (let i = 0; i < length; i++) {
let sum = 0;
for (let j = 0; j < numArrays; j++) { sum += newArrays[j][i] }
result[i] = sum / numArrays;
}
return result;
}; //вычисление среднего массива
function getCrossSpectrum(reX, imX, reY, imY) {
const n = reX.length;
const reXY = new Array(n); // действительная часть
const imXY = new Array(n); // мнимая часть
const power = new Array(n); // мощность кросс-спектра
let phase = new Array(n); // фаза (радианы)
for (let i = 0; i < n; i++) {
// Комплексное сопряжение спектра imY* (мнимая часть спектра силы сразу берём с обратным знаком)
const a = reX[i], b = imX[i]; // X(f) = a + ib
const c = reY[i], d = -imY[i]; // Y(f) = c + id
// Почастотное умножение
reXY[i] = a * c - b * d; // действительная часть
imXY[i] = b * c + a * d; // мнимания часть
// Мощность кросс-спектра
power[i] = reXY[i] * reXY[i] + imXY[i] * imXY[i];
// Фаза = arctan(Im/Re)
phase[i] = Math.atan2(imXY[i], reXY[i]);
}
return { reXY, imXY, power, phase };
}
function getFrequencyResponse(reX, imX, reY, imY) {
const n = reX.length;
const magnitude = new Array(n); // амплитуда (линейная)
const magnitudeDB = new Array(n); // амплитуда (дБ)
let phase = new Array(n); // фаза (радианы)
const RAD_TO_DEG = 180 / Math.PI; //перевод радиан в градусы
for (let i = 0; i < n; i++) {
const a = reX[i], b = imX[i]; // X(f) = a + ib
const c = reY[i], d = imY[i]; // Y(f) = c + id
// Знаменатель
const denominator = c * c + d * d; // знаменатель исходный вариант
// const denominator = a * a + b * b; // знаменатель вариант 2
// Избегаем деления на ноль
if (Math.abs(denominator) < 1e-10) {
magnitude[i] = 0;
magnitudeDB[i] = -Infinity;
phase[i] = 0;
continue;
}
// Почастотное умножение
const reH = (a * c + b * d) / denominator; // действительная часть
const imH = (c * b - a * d) / denominator; // мнимая часть исходный вариант
// const imH = (a * d - b * c) / denominator; // мнимая часть вариант 2 (поменял местами вычитаемые)
// Амплитуда
magnitude[i] = Math.sqrt(reH * reH + imH * imH);
// Амплитуда в дБ
magnitudeDB[i] = 20 * Math.log10(magnitude[i] / 1e-6);
// Приводим фазу сразу к диапазону [0°, 360°]
phase[i] = Math.atan2(imH, reH) * RAD_TO_DEG % 360;
if (phase[i] < 0) { phase[i] += 360 };
}
return { magnitude, magnitudeDB, phase };
}