Модель фильтра Калмана 2 порядка — различия между версиями

Материал из SRNS
Перейти к: навигация, поиск
(С моделью фазы)
 
(не показаны 16 промежуточных версий 1 участника)
Строка 14: Строка 14:
 
|repository    =  
 
|repository    =  
 
|category1    = Статистическая радиотехника
 
|category1    = Статистическая радиотехника
|category2    = Фазовые измерения
+
|category2    = Оценивание частоты
|category3    = Переходные процессы
+
|category3    = Оценивание задержки
 
|category4    =  
 
|category4    =  
|category5    =
+
|category5    =  
 
|category6    =   
 
|category6    =   
 
|category7    =   
 
|category7    =   
Строка 61: Строка 61:
 
for k = 1:K
 
for k = 1:K
 
     Ud = f(Xextr, Xist);  % Дискриминатор
 
     Ud = f(Xextr, Xist);  % Дискриминатор
     Sd = f(A_IQ);      % Критизна дискриминационной характеристики
+
     Sd = f(A_IQ);      % Крутизна дискриминационной характеристики
 
     Xest = Xextr + Ko*Ud/Sd; % Вектор оценок на k-й интервал
 
     Xest = Xextr + Ko*Ud/Sd; % Вектор оценок на k-й интервал
 
     Xextr = F*Xest;        % Экстраполяция на интервал k+1
 
     Xextr = F*Xest;        % Экстраполяция на интервал k+1
Строка 81: Строка 81:
 
qcno_ist = 45*ones(1,K); % // SNR
 
qcno_ist = 45*ones(1,K); % // SNR
  
H = 20; % Hz, полоса
+
H = 3; % Hz, полоса
  
 
Xextr = [0; 0; 0]; % Вектор экстраполяций
 
Xextr = [0; 0; 0]; % Вектор экстраполяций
Строка 119: Строка 119:
 
EpsTau = nan(1, K);
 
EpsTau = nan(1, K);
  
nI = stdn_IQ(k)*randn(1,K); % // I-comp noise
+
nI = stdn_IQ.*randn(1,K); % // I-comp noise
nQ = stdn_IQ(k)*randn(1,K); % // Q-comp noise
+
nQ = stdn_IQ.*randn(1,K); % // Q-comp noise
  
w = 0;
+
w = 0; Isum = 0; Qsum = 0; Iold = 1; Qold = 0;
 
for k = 1:K
 
for k = 1:K
  
Строка 129: Строка 129:
 
     EpsW(k) = Xist(2) - Xextr(2);
 
     EpsW(k) = Xist(2) - Xextr(2);
 
     EpsTau(k) = 0;
 
     EpsTau(k) = 0;
     [A_IQ(k) qcno] = qcno_change(qcno_ist(k), stdn_IQ(k), Tc);
+
 
 +
     qcno = 10.^(qcno_ist(k)/10);
 +
    A_IQ(k) = stdn_IQ(k) .* sqrt(2 * qcno * Tc);
 +
 
 
     A_IQ_eff(k) = A_IQ(k)*sinc(EpsW(k)*Tc/2 /pi)*ro(EpsTau(k));
 
     A_IQ_eff(k) = A_IQ(k)*sinc(EpsW(k)*Tc/2 /pi)*ro(EpsTau(k));
 
     mI = A_IQ_eff(k) * cos(EpsW(k)*Tc/2 + EpsPhi(k));
 
     mI = A_IQ_eff(k) * cos(EpsW(k)*Tc/2 + EpsPhi(k));
Строка 135: Строка 138:
 
     I(k) = mI + nI(k);
 
     I(k) = mI + nI(k);
 
     Q(k) = mQ + nQ(k);
 
     Q(k) = mQ + nQ(k);
      
+
     Isum = Isum + I(k);
 +
    Qsum = Qsum + Q(k);
 +
 
 
     Xextr = Fincorr * Xextr; % Набег фазы в корреляторе к концу накопления     
 
     Xextr = Fincorr * Xextr; % Набег фазы в корреляторе к концу накопления     
 
      
 
      
 
     w = w + 1;
 
     w = w + 1;
     if w == (Tf/Tc)                   
+
     if w == fix(Tf/Tc)                   
         Ud = f(I(k), Q(k), I(k-1), Q(k-1), ...); % Дискриминатор, работающий по некоторому набору квадратур
+
         %    Ud = f(I(k), Q(k), I(k-1), Q(k-1), ...);     % Дискриминатор
         Sd = f(A_IQ);                           % Крутизна дискриминационной характеристики
+
         Ud = (I(k)*Qold - Q(k)*Iold);
 +
        %    Sd = f(A_IQ);             % Крутизна дискриминационной характеристики
 +
        Sd = Tc*(A_IQ(k)*Tf/Tc)^2 * 1.3;
 
         Xest = Xextr + Ko*Ud/Sd;                % Вектор оценок на очередной интервал фильтра
 
         Xest = Xextr + Ko*Ud/Sd;                % Вектор оценок на очередной интервал фильтра
 
         Xextr = Ff*Xest;                        % Экстраполяция на следующий интервал
 
         Xextr = Ff*Xest;                        % Экстраполяция на следующий интервал
 
         w = 0;
 
         w = 0;
 +
        Iold = Isum; Isum = 0;
 +
        Qold = Qsum; Qsum = 0;
 
     end
 
     end
  
     Xist = Fc*Xist + [0; 1]*nIst(k)*stdIst; % Здесь может быть любая другая модель изменения истинного вектора состояния
+
     Xist = Fc*Xist + [0; 0; 1]*nIst(k)*stdIst; % Здесь может быть любая другая модель изменения истинного вектора состояния
     Iold = I(k);
+
     if ~mod(k,fix(K/10))
    Qold = Q(k);
+
        fprintf('Progress: %.0f%%\n', 100*k/K);
 +
    end
 
end
 
end
  
Строка 155: Строка 165:
  
 
== Результаты моделирования при параметрах по-умолчанию ==
 
== Результаты моделирования при параметрах по-умолчанию ==
 +
 +
<gallery>
 +
Image:20111012_FLL_1.png|Ошибка слежения за частотой для ЧАП из примера
 +
</gallery>
  
 
== См. также ==
 
== См. также ==

Текущая версия на 17:05, 22 июля 2016

Список всех моделей
Crystal Clear action run.png
Модель фильтра Калмана 2 порядка
Описание Модель фильтра Калмана 2 порядка на примере ФАП
Автор(ы) Korogodin (Korogodinобсуждение)
Последняя версия 1.0 (12.10.2011)
Загрузить no link
Хранилище no link
Категории Статистическая радиотехника, Оценивание частоты, Оценивание задержки


Содержание

[править] Описание модели

Модель фильтра Калмана 2 порядка, например, используемого в ЧАП. В данный момент приведен листинг только для коэффициентов установившегося режима. Следует привести пример с уравнениями Рикатти.

[править] Листинг

Ниже приведен листинг при использовании коэффициентов установившегося режима. Изложение следует дополнить уравнениями Рикатти - для честного соответствия заголовку.

Коэффициенты передачи непрерывных следящих систем второго и третьего порядка

[править] Без моделирования фазы (первообразной от главного оцениваемого параметра)

Tmod = 300; % Время моделирования
Tc = 0.001; % Период работы фильтров
K = fix(Tmod/Tc);

Xextr = [0; 0]; % Вектор экстраполяций

F = [1 Tc
     0 1  ]; % Переходная матрица

H = 20; % Hz, полоса

Ko = nan(2,1); % Вектор-столбец коэффициентов фильтра
Ko(2) = 2*16/9*H^2; % Коэффициенты непрерывной системы в установившемся режиме
Ko(1) = sqrt(2*Ko(2));

Ko = Ko*Tc; % Переход к коэффициентам дискретной системы

Xist = [0; 0]; % Истинный вектор состояния
stdIst = 10; nIst = randn(1,K);
for k = 1:K
    Ud = f(Xextr, Xist);  % Дискриминатор
    Sd = f(A_IQ);      % Крутизна дискриминационной характеристики
    Xest = Xextr + Ko*Ud/Sd; % Вектор оценок на k-й интервал
    Xextr = F*Xest;         % Экстраполяция на интервал k+1

    Xist = F*Xist + [0; 1]*nIst(k)*stdIst; % Здесь может быть любая другая модель изменения истинного вектора состояния
end

[править] С моделью фазы

Для моделирования ЧАП с статистическими эквивалентами корреляторов необходимо учитывать разность набегов фазы приходящего сигнала и фазы опорного колебания в корреляторе. Тогда удобно использовать следующую модель:

Tmod = 300; % Время моделирования
Tf = 0.005; % Период работы фильтров
Tc = 0.001; % Период интегрирования в корреляторе
K = fix(Tmod/Tc);

qcno_ist = 45*ones(1,K); % // SNR

H = 3; % Hz, полоса

Xextr = [0; 0; 0]; % Вектор экстраполяций

Ff = [1 0  0
      0 1  Tf
      0 0  1]; % Переходная матрица для модели частоты (в темпе фильтра)
 
Fc = [1 Tc Tc^2/2
      0 1  Tc
      0 0  1]; % Переходная матрица для модели частоты (в темпе коррелятора)
 
Fincorr =  [1 Tc 0
            0 1  0
            0 0  1]; % Переходная матрица набега фазы в корреляторе

Ko = nan(3,1); % Вектор-столбец коэффициентов фильтра
Ko(3) = 2*16/9*H^2; % Коэффициенты непрерывной системы в установившемся режиме
Ko(2) = sqrt(2*Ko(3));
Ko(1) = 0;    % Это не баг, это фича! Из-за этого нуля система, на самом деле, - второго порядка

Ko = Ko*Tf; % Переход к коэффициентам дискретной системы

Xist = [0; 0; 0]; % Истинный вектор состояния
stdIst = 10; nIst = randn(1,K);

stdn_IQ = ones(1,K)*8; % СКО шума квадратурных сумм

A_IQ = nan(1,K); % // Memory allocation
A_IQ_eff = nan(1,K);

I = nan(1,K); % // Memory allocation
Q = nan(1,K);

EpsPhi = nan(1, K); % // Memory allocation
EpsW = nan(1, K);
EpsTau = nan(1, K);

nI = stdn_IQ.*randn(1,K); % // I-comp noise
nQ = stdn_IQ.*randn(1,K); % // Q-comp noise

w = 0; Isum = 0; Qsum = 0; Iold = 1; Qold = 0;
for k = 1:K

    % // Расчет стат.эквивалентов корреляционных сумм
    EpsPhi(k) = mod(Xist(1) - Xextr(1),2*pi);
    EpsW(k) = Xist(2) - Xextr(2);
    EpsTau(k) = 0;

    qcno = 10.^(qcno_ist(k)/10);
    A_IQ(k) = stdn_IQ(k) .* sqrt(2 * qcno * Tc);

    A_IQ_eff(k) = A_IQ(k)*sinc(EpsW(k)*Tc/2 /pi)*ro(EpsTau(k));
    mI = A_IQ_eff(k) * cos(EpsW(k)*Tc/2 + EpsPhi(k));
    mQ = - A_IQ_eff(k) * sin(EpsW(k)*Tc/2 + EpsPhi(k));
    I(k) = mI + nI(k);
    Q(k) = mQ + nQ(k);
    Isum = Isum + I(k);
    Qsum = Qsum + Q(k);

    Xextr = Fincorr * Xextr; % Набег фазы в корреляторе к концу накопления    
   
    w = w + 1;
    if w == fix(Tf/Tc)                  
        %     Ud = f(I(k), Q(k), I(k-1), Q(k-1), ...);      % Дискриминатор
        Ud = (I(k)*Qold - Q(k)*Iold);
        %     Sd = f(A_IQ);             % Крутизна дискриминационной характеристики
        Sd = Tc*(A_IQ(k)*Tf/Tc)^2 * 1.3;
        Xest = Xextr + Ko*Ud/Sd;                 % Вектор оценок на очередной интервал фильтра
        Xextr = Ff*Xest;                         % Экстраполяция на следующий интервал
        w = 0;
        Iold = Isum; Isum = 0;
        Qold = Qsum; Qsum = 0;
    end

    Xist = Fc*Xist + [0; 0; 1]*nIst(k)*stdIst; % Здесь может быть любая другая модель изменения истинного вектора состояния
    if ~mod(k,fix(K/10))
        fprintf('Progress: %.0f%%\n', 100*k/K);
    end
end

[править] Результаты моделирования при параметрах по-умолчанию

[править] См. также

Персональные инструменты
Пространства имён

Варианты
Действия
SRNS Wiki
Рабочие журналы
Приватный файлсервер
QNAP Сервер
Инструменты