Оцените передаточную функцию неизвестной системы

Можно оценить передаточную функцию неизвестной системы на основе измеренных входных и выходных данных системы.

В DSP System Toolbox™ можно оценить передаточную функцию системы с помощью dsp.TransferFunctionEstimator Система object™ в MATLAB^® и блок Discrete Transfer Function Estimator в Simulink^®. Отношение между входом x и выходом y моделируется линейной, независимой от времени передаточной функцией _Txy. Передаточная функция является отношением взаимной спектральной плотности мощности x и y, _Pyx, к спектральной плотности мощности x, _Pxx:

$T_{x y} (f) = \frac{P_{y x} (f)}{P_{x x} (f)}$

dsp.TransferFunctionEstimator объект и блок Discrete Transfer Function Estimator используют усредненный метод периодограммы валлийцев, чтобы вычислить _Pxx и _Pxy. Для получения дополнительной информации об этом методе смотрите Спектральный анализ.

Когерентность

Когерентность или квадрат когерентности, между x и y задана как:

$C_{x y} (f) = \frac{{| P_{x y} |}^{2}}{P_{x x} * P_{y y}}$

Функция когерентности оценивает степень, до которой можно предсказать y от x. Значение когерентности находится в области значений 0 ≤ _Cxy (f) ≤ 1. Если _Cxy = 0, вход x и выход y не связан. Значение _Cxy, больше, чем 0 и меньше чем 1 указывает на одно из следующего:

Измерения являются шумными.
Система нелинейна.
Выход y является функцией x и других входных параметров.

Когерентность линейной системы представляет дробную часть энергии выходного сигнала, которая произведена входом на той частоте. Для особой частоты 1 – _Cxy является оценкой дробной степени выхода, что вход не способствует.

Когда вы устанавливаете OutputCoherence свойство dsp.TransferFunctionEstimator к true, объект вычисляет выходную когерентность. В блоке Discrete Transfer Function Estimator, чтобы вычислить спектр когерентности, устанавливают флажок Output magnitude squared coherence estimate.

Оцените передаточную функцию в MATLAB

Открыть скрипт

Чтобы оценить передаточную функцию системы в MATLAB™, используйте dsp.TransferFunctionEstimator Система object™. Объект реализует средний модифицированный метод периодограммы валлийцев и использует измеренные входные и выходные данные для оценки.

Инициализируйте систему

Система является каскадом двух этапов фильтра: dsp.LowpassFilter и параллельная связь dsp.AllpassFilter и dsp.AllpoleFilter.

allpole = dsp.AllpoleFilter;
allpass = dsp.AllpassFilter;
lpfilter = dsp.LowpassFilter;

Задайте источник сигнала

Вход к системе является синусоидой с частотой 100 Гц. Частота дискретизации составляет 44,1 кГц.

sine = dsp.SineWave('Frequency',100,'SampleRate',44100,...
    'SamplesPerFrame',1024);

Создайте средство оценки передаточной функции

Чтобы оценить передаточную функцию системы, создайте dsp.TransferFunctionEstimator Системный объект.

tfe  = dsp.TransferFunctionEstimator('FrequencyRange','onesided',...
    'OutputCoherence', true);

Создайте график массивов

Инициализируйте два dsp.ArrayPlot объекты: один, чтобы отобразить ответ величины системы и другого, чтобы отобразить оценку когерентности между входом и выходом.

tfeplotter = dsp.ArrayPlot('PlotType','Line',...
    'XLabel','Frequency (Hz)',...
    'YLabel','Magnitude Response (dB)',...
    'YLimits',[-120 20],...
    'XOffset',0,...
    'XLabel','Frequency (Hz)',...
    'Title','System Transfer Function',...
    'SampleIncrement',44100/1024);
coherenceplotter = dsp.ArrayPlot('PlotType','Line',...
    'YLimits',[0 1.2],...
    'YLabel','Coherence',...
    'XOffset',0,...
    'XLabel','Frequency (Hz)',...
    'Title','Coherence Estimate',...
    'SampleIncrement',44100/1024);

По умолчанию ось X графика массивов находится в выборках. Чтобы преобразовать эту ось в частоту, установите свойство 'SampleIncrement' dsp.ArrayPlot возразите против Фс/1024. В этом примере это значение является 44100/1024, или 43.0664. Для двухстороннего спектра, XOffset свойство dsp.ArrayPlot объект должен быть [-Fs/2]. Частота варьируется по области значений [-Fs/2 Фс/2]. В этом примере график массивов показывает односторонний спектр. Следовательно, установите XOffset к 0. Частота варьируется по области значений [0 Фс/2].

Оцените передаточную функцию

Средство оценки передаточной функции принимает два сигнала: введите к фильтру 2D этапа и выходу фильтра 2D этапа. Вход к фильтру является синусоидой, содержащей аддитивный белый Гауссов шум. Шум имеет среднее значение нуля и стандартное отклонение 0,1. Средство оценки оценивает передаточную функцию фильтра 2D этапа. Выход средства оценки является частотной характеристикой фильтра, который является комплексным. Чтобы извлечь фрагмент величины этой комплексной оценки, используйте функцию abs. Чтобы преобразовать результат в дБ, примените коэффициент преобразования 20*log10 (величина).

for Iter = 1:1000
    input = sine() + .1*randn(1024,1);
    lpfout = lpfilter(input);
    allpoleout = allpole(lpfout);
    allpassout = allpass(lpfout);
    output = allpoleout + allpassout;
    [tfeoutput,outputcoh] = tfe(input,output);
    tfeplotter(20*log10(abs(tfeoutput)));
    coherenceplotter(outputcoh);
end

Первый график показывает ответ величины системы. Второй график показывает оценку когерентности между вводом и выводом системы. Когерентность в графике варьируется по области значений [0 1] как ожидалось.

Ответ величины Фильтра Используя fvtool

Фильтр является каскадом двух этапов фильтра - dsp.LowpassFilter и параллельная связь dsp.AllpassFilter и dsp.AllpoleFilter. Все объекты фильтра используются в их состоянии по умолчанию. Используя коэффициенты фильтра, выведите системную передаточную функцию и постройте частотную характеристику с помощью freqz. Ниже коэффициенты в [Цифре] [Логово] формат:

Весь фильтр полюса - [1 0] [1 0.1]
Весь фильтр передачи - [0.5 - 1/sqrt (2) 1] [1 - 1/sqrt (2) 0.5]
Фильтр lowpass - Определяет коэффициенты с помощью следующих команд:

lpf = dsp.LowpassFilter;
Coefficients = coeffs(lpf);

Coefficients.Numerator дает коэффициенты в формате массивов. Математическую деривацию полной системной передаточной функции не показывают здесь. Если вы выводите передаточную функцию, запускаете fvtool, и вы видите частотную характеристику ниже:

Ответ величины, что fvtool показывает соответствиям ответ величины что dsp.TransferFunctionEstimator объектные оценки.

Оцените передаточную функцию в Simulink

Чтобы оценить передаточную функцию системы в Simulink, используйте блок Discrete Transfer Function Estimator. Блок реализует средний модифицированный метод периодограммы валлийцев и использует измеренные входные и выходные данные для оценки.

Система является каскадом двух этапов фильтра: фильтр lowpass и параллельная связь фильтра allpole и фильтра allpass. Вход к системе является синусоидой, содержащей аддитивный белый Гауссов шум. Шум имеет среднее значение нуля и стандартное отклонение 0,1. Вход к средству оценки является системным входом и системой выход. Выход средства оценки является частотной характеристикой системы, которая является комплексной. Чтобы извлечь фрагмент величины этой комплексной оценки, используйте блок Abs. Чтобы преобразовать результат в дБ, система использует блок dB (1 ohm).

Snapshot of the transfer function estimator model. Source blocks on the left followed by the filter cascade blocks, transfer function estimator, Abs block, dB Conversion block, and two Array Plot blocks.

Откройте и смотрите модель

Чтобы открыть модель, введите ex_transfer_function_estimator в командной строке MATLAB.

Вот настройки блоков в модели.

Блок	Изменения параметра	Цель блока
Sine Wave	Sample time к `1/44100` Samples per frame к `1024`	Сигнал синусоиды с частотой на уровне 100 Гц
Random Source	Source type к `Gaussian` Variance к 0,01 Sample time к 1/44100 Samples per frame к 1 024	Блок Random Source генерирует случайный шумовой сигнал со свойствами, заданными через диалоговое окно блока
Lowpass Filter	Никакое изменение	Фильтр lowpass
Allpole Filter	Никакое изменение	Allpole фильтруют с коэффициентами `[1 0.1]`
Discrete Filter	Numerator к `[0.5 -1/sqrt(2) 1]` Denominator к `[1 -1/sqrt(2) 0.5]`	Allpass фильтруют с коэффициентами `[-1/sqrt(2) 0.5]`
Discrete Transfer Function Estimator	Frequency range к `One-sided` Number of spectral averages к 8	Средство оценки передаточной функции
Abs	Никакое изменение	Извлекает информацию о величине из выхода средства оценки передаточной функции
Первый блок Array Plot	Нажмите View: Выберите Style и установите Plot type на `Line`. Выберите Configuration Properties: От вкладки Main, набор Sample increment к `44100/1024` и X-offset к `0`. Во вкладке Display задайте Title как `Magnitude Response of the System in dB`, X-label как `Frequency (Hz)`, и Y-label как `Amplitude (dB)`	Показывает ответ величины системы
Второй блок Array Plot	Нажмите View: Выберите Style и установите Plot type на `Line`. Выберите Configuration Properties: От вкладки Main, набор Sample increment к `44100/1024` и X-offset к `0`. Во вкладке Display задайте Title как `Coherence Estimate`, X-label как `Frequency (Hz)`, и Y-label как `Amplitude`	Показывает оценку когерентности

По умолчанию x - ось графика массивов находится в выборках. Чтобы преобразовать эту ось в частоту, параметр Sample increment устанавливается на Fs/1024. В этом примере этим значением является 44100/1024, или 43.0664. Для двухстороннего спектра параметром X-offset должен быть –Fs/2. Частота варьируется по области значений [-Fs/2 Fs/2]. В этом примере график массивов показывает односторонний спектр. Следовательно, X-offset установлен в 0. Частота варьируется по области значений [0 Fs/2].

Запустите модель

Output of the two Array Plot blocks are placed next to each other. First Array Plot block on the left shows the magnitude response of the system in dB. The second Array Plot block on the right shows the coherence estimate. The x-axis for both plots shows the frequencies in Hz, ranging from 0 to 22500 Hz.

Документация