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

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

В 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)}$

The dsp.TransferFunctionEstimator объект и Discrete Transfer Function Estimator блок используют усредненный метод периодограммы Welch для вычисления _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™. Объект реализует среднее значение модифицированный метод периодограммы Welch и использует измеренные входные и выходные данные для оценки.

Инициализация системы

Система представляет собой каскад из двух ступеней фильтра: 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 объект к Fs/1024. В этом примере это значение является 44100/1024 или 43.0664. Для двустороннего спектра, XOffset свойство dsp.ArrayPlot объект должен быть [-Fs/2]. Частота изменяется в области значений [-Fs/2 Fs/2]. В этом примере график массива показывает односторонний спектр. Следовательно, установите XOffset в 0. Частота изменяется в области значений [0 Fs/2].

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

Блок оценки передаточной функции принимает два сигнала: вход в двухкаскадный фильтр и выход двухкаскадного фильтра. Вход фильтра представляет собой синусоиду, содержащую аддитивный белый Гауссов шум. Шум имеет среднее значение нуля и стандартное отклонение 0,1. Оценщик оценивает передаточную функцию двухэтапного фильтра. Выходом оценщика является частотная характеристика фильтра, который является комплексным. Чтобы извлечь фрагмент величины этой комплексной оценки, используйте функцию 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. Ниже приведены коэффициенты в формате [Num] [Den]:

Весь шестовый фильтр - [1 0] [1 0,1]
Весь пропускной фильтр - [0,5 -1/sqrt (2) 1] [1 -1/sqrt (2) 0,5]
Lowpass фильтр - Определите коэффициенты с помощью следующих команд:

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

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

Величина ответ, который fvtool показов, совпадает с величиной ответом, который dsp.TransferFunctionEstimator оценки объектов.

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

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

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

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

Чтобы открыть модель, введите 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 - 1024	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].

Запуск модели

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