Цель работы: анализ динамики частотных характеристик физической системы и особенностей представления работы сисстемы в ЭВМ.
Исследование имеет целью разработку алгоритма, который в реальном времени по данным измерения дискретизованного с заданной частотой выходного сигнала позволяет определить номер режима в котором работает система.
Работа системы в стационарном режиме имитируется программой dtmp.m. Номера режимов и соответствующие значения параметров системы, при которых надлежит провести исследование составляют спецификацию индивидуального задания.
Инстументарий для выполнения анализа представлен в виде программ, написанных на языке Matlab 3.5.
Итогом исследования является вариант параметрического описания набора цифровых фильтров, применение которых к сигналам системы позволит определить состояние системы, и алгоритм принятия решения по результатам фильтрации.
| Номер варианта |
sens |
tt мксек |
pp мксек |
Частота дискр-ии |
Н-ра режимов работы |
| 1 |
0.04 |
30 |
120 |
8000 |
0,..,3 |
| 2 |
0.2 |
40 |
140 |
16000 |
3,..,6 |
| 3 |
0.02 |
20 |
60 |
22400 |
1,..,4 |
| 4 |
0.03 |
25 |
40 |
8000 |
6,..,9 |
| 5 |
0.04 |
15 |
90 |
16000 |
0,..,3 |
| 6 |
0.06 |
60 |
180 |
22400 |
3,..,6 |
| 7 |
0.08 |
15 |
50 |
8000 |
1,..,4 |
| 8 |
0.01 |
10 |
60 |
16000 |
6,..,9 |
| 9 |
0.1 |
25 |
100 |
22400 |
0,..,3 |
| 10 |
0.01 |
10 |
70 |
8000 |
3,..,6 |
| 11 |
0.04 |
30 |
120 |
16000 |
1,..,4 |
| 12 |
0.2 |
40 |
140 |
22400 |
6,..,9 |
| 13 |
0.02 |
20 |
60 |
8000 |
0,..,3 |
| 14 |
0.03 |
25 |
40 |
16000 |
3,..,6 |
| 15 |
0.04 |
15 |
90 |
22400 |
1,..,4 |
| 16 |
0.06 |
60 |
180 |
8000 |
6,..,9 |
| 17 |
0.04 |
30 |
120 |
16000 |
0,..,3 |
| 18 |
0.01 |
10 |
60 |
16000 |
6,..,9 |
| 19 |
0.1 |
25 |
100 |
22400 |
0,..,3 |
| 20 |
0.01 |
10 |
70 |
8000 |
3,..,6 |
| 21 |
0.04 |
30 |
120 |
16000 |
1,..,4 |
| 22 |
0.2 |
40 |
140 |
22400 |
6,..,9 |
| 23 |
0.02 |
20 |
60 |
8000 |
0,..,3 |
| 24 |
0.03 |
25 |
40 |
16000 |
3,..,6 |
| 25 |
0.04 |
15 |
90 |
22400 |
1,..,4 |
function cont=dtmp(t)
% global tt pp
% floor(a)) - округляет элементы a до
% ближайшего в направлении к (-inf) значения
t1=t-(tt+pp)*floor(t./(tt+pp));
% t1 - остаток от деления t на период (tt+pp)
cont=(t1<=tt).*(0.9+t1/tt*0.1);
% cont = 0, если доля t1 в (tt+pp) больше tt,
% иначе - cont растет от 0.9 до 1
cont=cont(:);
**************************************************************
Примеры исследования спектра и верификации предлагаемых алгоритмов
приведены в виде последовательности скриптов
tst1
tst2
filt1
filt2
filttst
Последовательность исполнения существенна.
В отчет включить параметры фильтров (тип,частотный диапазон, ко-
эффициенты в формате long) и константы порогов принятия решений
*************************************************************
tst1.m
*************************************************************
% tst1 - моделирует работу стенда с аналоговым
% представлением данных - на экране осцилографа
% dtmp - функция, моделирующая
% температурно-зависимый генератор импульсов
clear all % Очистка памяти
diary tst1.txt % Включение записи
% текстовых обменов в протокольный файл tst1.txt
global tt pp sens % Описание глобальных переменных, используемых
% в tst1 и tst2
sens=0.02; % Условная характеристика термочувствительности
% генератора, определяемая типом транзистора
% и режимом его работы
Tdiap=2 % При температурных градациях i=0:(Tdiap-1)
% исследуем спектр колебаний
N=256; % Длина блока, допустимого по ограничениям памяти
tt=35*10^(-6); % Длина токовой фазы
pp0=3*tt; % Длина бестоковой фазы при нулевой температуре
fd=10^6; % Частота дискретизации для имитации осцилограммы
% fd=8000; % Частота дискретизации для зв. карты
% fd=16000; % Частота дискретизации для зв. карты
fs=fd/2 % Частота Найквиста
step=1/fd; % Интервал дискретизации
t=0:step:step*(N-1); % Массив временных отсчетов для осцилограммы
at=length(t); % Его длина
yy=[]; % Инициализация массива данных амплитуд
for ts=0:(Tdiap-1) % Начало цикла по температуре
pp=(1-sens*ts)*pp0 % Изменяем настройку:
% в соответствии с ростом температуры
% уменьшается длина бестоковой фазы
y=dtmp(t+rand); % Моделируем блок данных от генератора со
% случайной фазой
%y=detrend(y,0);
yy=[yy,y]; % Добавляем его в массив
end % Конец цикла по температуре
xp=max(t);yp=max(yy(:,1)*7/8);
% plot(t(1:N/8),yy(1:N/8,:)),title(' График во временной области ')
plot(t,yy),title(' График во временной области ')
text(xp*0.6,yp,['Чувствительность =',num2str(sens)]),grid,pause
% Графики во временной области
% частоты
save tst1
*************************************************************
*************************************************************
tst2.m
*************************************************************
% tst2 для dtmp - Функция, моделирующая
% температурно-зависимый генератор импульсов
clear all % Очистка памяти
load tst1 % Загрузка глобальных переменных
diary tst2.txt % Включение записи
% текстовых обменов в протокольный файл tst2.txt
global tt pp sens % глобальные переменные
% sens - Условная характеристика термочувствительности
% генератора, определяемая типом транзистора и режимом его работы
Tdiap=2 % - при температурных градациях i=0:(Tdiap-1)
% исследуем спектр колебаний
% N=256; длина блока, допустимого по ограничениям памяти
fd=8000; % Частота дискретизации для зв. карты - Версия 1
% fd=16000; % Частота дискретизации для зв. карты - Версия 2
fs=fd/2 % Частота Найквиста
step=1/fd; % Интервал дискретизации
t=0:step:(N-1)*step; % Массив временных отсчетов
at=length(t); % Его длина
gg=floor(log(at)/log(2));
% gg - наименьший показатель степени, для которой
% 2^gg укладывается на считанном блоке отсчетов
ffl=2^(gg+1); % Длина FFT, возможно, с дополнительными нулями
rr=1:ffl/2; % Берем значимую область частот,
% удаляя симметричные данные
rr=rr(:); % Преобразуем в вектор-колонку
lr=length(rr); % Длина этого вектора
axe=2:lr; % Это нужно для удаления нулевой компоненты,
% которая сответствует 1-му члену ряда Фурье
aftt=[]; % Инициализация массива мощностей частотных компонент
q12=[]; % Инициализация массива мощностей и
% номеров частотных компонент
yy=[]; % Инициализация массива данных амплитуд
for ts=0:(Tdiap-1) % Начало цикла по температуре
pp=(1-sens*ts)*pp0 % Изменяем настройку:
% в соответствии с ростом температуры
% уменьшается длина бестоковой фазы
y=dtmp(t+rand); % Моделируем блок данных от генератора
% y=detrend(y,0); % Эта функция удалит постоянную составляющую
yy=[yy,y]; % Добавляем блок данных в массив
ft=fft(y,ffl); % Находим Фурье-представление данных
rr=ft(axe); % Удалим нулевую компоненту спектра
pow=rr.*conj(rr); % Находим мощности спeктральных компонент
aftt=[aftt,pow]; % Добавляем вектор мощностей спектра в массив afft
% для построения графика
[q1,q2]=sort(pow); % Сортируем pow по возрастанию мощностей,
% сохраняя номера в q2
q12=[q12,q1,q2]; % Добавим в q12 номера и энергии ведущих компонент
end % Конец цикла по температуре
nmax=q12(lr-1:-1:lr-6,:); % Массив номеров и энергий
% максимальных спектральных компонент
xp=max(t); yp=max(yy(:,1)*7/8);
% plot(t(1:N/3),yy(1:N/3,:)),title(' График во временной области ')
% Этот график использует только первую треть данных
plot(t(1:N),yy(1:N,:)),title(' График во временной области ')
text(xp*0.6,yp,['Чувствительность =',num2str(sens)]),grid,pause
% Графики во временной области
% но на частичном временном отрезке
plot(axe*fd/ffl,aftt),title('Спектральные плотности '),grid,pause
% Графики во частотной области
disp(['Цена деления по шкале частот=', num2str(fd/lr)]);
disp(['Частота Найквиста=', num2str(fs)]);
disp(['Мощности и номера ведущих спектральных компонент']);
nmax, pause
disp(['Ведущие спектральные частоты (гц) ']);
fmax=fd/ffl*nmax(:,2:2:2*Tdiap) % Переводим номера ведущих
% компонент в соответствующие частоты
save tst2
*************************************************************
*************************************************************
filt1.m
*************************************************************
load tst2 % Загрузка переменных, вычисленных в tst2
diary filt1.txt % Включение записи
% текстовых обменов в протокольный файл filt1.txt
global tt pp sens % глобальные переменные
% Спроектируем полосовай фильтр Баттеруорта
% пропускающий частоту, которая доминирует на
% первом температурном уровне и не проявляется на
% других
% Данные о доминирующих частотах
% 1.0e+003 *
% 0.8594 0.7500
% 0.8438 0.0156
% 0.0156 1.5000
% 1.7188 0.7344
% 0.8750 0.7656
% 1.7031 1.4844
format long
[b1,a1]=butter(5,[800,900]/4000)
format
% Проверим, что фильтр действительно пропускает
% значительно лучше первую кривую, чем вторую
u=[];
for ts=1:2 % Начало цикла по температуре
uu=yy(:,ts);
uu=filter(b1,a1,uu);
u=[u,uu];
plot(t(1:N),u(1:N,:)),title(' График во временной области '),grid,pause
end
save tst2
% Теперь можно дать колличественную оценку результату фильтрации
id=mean(abs(u))
% Напоминание: использованы функции MEAN - "среднее"
$ и ABS - "абс. значение".
% MEAN Average or mean value. For vectors, MEAN(X) is the mean
% value of the elements in vector X. For matrices, MEAN(x)
% is a row vector containing the mean value of each column.
% ABS ABS(X) is the absolute value of the elements of X. When
% X is complex, ABS(X) is the complex modulus (magnitude) of
% the elements of X. See also ANGLE.
*************************************************************
*************************************************************
filt2.m
*************************************************************
load tst2 % Загрузка переменных, вычисленных в tst2
diary filt1.txt % Включение записи
% текстовых обменов в протокольный файл filt1.txt
global tt pp sens % глобальные переменные
% Спроектируем фильтр на частоту, доминирующую на
% втором температурном уровне и не проявляющуюся
% на других
% Данные о доминирующих частотах
% 1.0e+003 *
% 0.8594 0.7500
% 0.8438 0.0156
% 0.0156 1.5000
% 1.7188 0.7344
% 0.8750 0.7656
% 1.7031 1.4844
format long
[b2,a2]=butter(5,[1450,1550]/4000)
format
u=[];
for ts=1:2 % Начало цикла по температуре
uu=yy(:,ts);
uu=filter(b2,a2,uu);
u=[u,uu];
plot(t(1:N),u(1:N,:)),title(' График во временной области '),grid,pause
end
save tst2
% Теперь можно дать колличественную оценку результату фильтрации
id=mean(abs(u))
*************************************************************
*************************************************************
filttst.m
*************************************************************
load tst2 % Загрузка переменных, сохраненных
% после вычислений в filt1 и filt2 в tst2
diary filttst.txt % Включение записи
% текстовых обменов в протокольный файл
global tt pp sens % глобальные переменные
% Теперь можно проверить алгоритм
% распознавания по результату фильтрации
for ts=0:(Tdiap-1) % Начало цикла по температуре
pp=(1-sens*ts)*pp0 % Изменяем настройку:
% в соответствии с ростом температуры
% уменьшается длина бестоковой фазы
y=dtmp(t+rand); % Моделируем блок данных от генератора
% y=detrend(y,0); % Иногда это улучшает работу
yy=[yy,y]; % Добавляем его в массив
end % Конец цикла по температуре
u1=[];
% Оба фильтра применим к отсчетам 1-го режима
uu=yy(:,1); % Считаем данные 1-го режима
u1=[u1,filter(b1,a1,uu)]; % Пропустим через фильтры
u1=[u1,filter(b2,a2,uu)];
plot(t(1:N),u1(1:N,:)),title(' График во временной области: 1-й режим'),grid,pause
u2=[]; % Оба фильтра применим к отсчетам 2-го режима
uu=yy(:,2); % Считаем данные 2-го режима
u2=[u2,filter(b1,a1,uu)]; % Пропустим через фильтры
u2=[u2,filter(b2,a2,uu)];
plot(t(1:N),u2(1:N,:)),title(' График во временной области: 2-й режим'),grid,pause
save filttst
% Теперь можно дать колличественную оценку результату фильтрации
id1=mean(abs(u1))
id2=mean(abs(u2))
*************************************************************
Теперь в качестве порога для принятия решения в пользу одного из решений достаточно провести раздел по значению bound=(id1+id2)/2.