1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

%Шаг 1. Генерируем шкалу времени от 0 до 1 секунды с шагом в 1/16384 секунды
(получится всего 16385 шагов). Массив t из 16385 элементов содержит значения времени.
t = 0:1/16384:1;

%Шаг 2. Генерим массивы амплитудами сигналов. В результате получим массивы
U25..U5555m12, в каждом из которых 16385 значения сигнала.
U25      = 1   *sin(2*pi *25.2    *t);
U75      = 1   *sin(2*pi* 75.5    *t);
U325     = 1   *sin(2*pi* 325.3   *t);
U420     = 0.1 *sin(2*pi* 420     *t);
U420p8   = 0.07*sin(2*pi*(420+8)  *t);
U420m8   = 0.07*sin(2*pi*(420-8)  *t);
U780     = 0.1 *sin(2*pi* 780     *t);
U780p12  = 0.07*sin(2*pi*(780+12) *t);
U780m12  = 0.07*sin(2*pi*(780-12) *t);
U5555    = 0.1 *sin(2*pi* 5555    *t);
U5555p12 = 0.07*sin(2*pi*(5555+12)*t);
U5555m12 = 0.07*sin(2*pi*(5555-12)*t);

%Шаг 3. Суммируем все сигналы между собой. В результате получим массив 
x из 16385 значений.
x = U25+U75+U325+U420+U420p8+U420m8+U780+U780p12+U780m12+U5555+U5555p12+U5555m12;

%Шаг 4. Для того, чтобы сигнал отображался "красиво", вводим оконную функцию
Хемминга. В результате получим массив оконного распределения w из 16385 значений.
Операция ".'" нужна для транспонирования оконного массива (он был столбцом, нужно
получить строку).
w = hamming(16385).';

%Шаг 5. Чтобы было "совсем реально", добавим немного шума. Здесь "size(t)" соответствует
"1, 16385", т. е. вызов randn(size(t)) это то же самое, что вызов randn(1, 16385).
x = x + 0.001*randn(size(t));

%Шаг 6. Применим оконную функцию к сигналу x, в результате получим массив z из
16385 значений. Здесь операция ".*" означает "перемножить каждый элемент массива
w с соответствующим элементом массива x"
z = w.*x;

%Шаг 7. Делаем дискретное преобразование фурье (ДПФ) от сигнала z по 16384 точкам.
В результате получим массив комплексных чисел Y (всего 16384 комплексных числа)
Y = fft(z,16384);

%Шаг 8. Вычисляем спектральную плотность по мощности, короче говоря - спектр
(чего мы и добивались). В результате получим массив Pyy из 16384 значений. Здесь
коэффициенты 5000 и 10000 подобраны таким образом, чтобы сместить значение
десятичного логарифма в диапазон значений целого беззнакового 16-битного числа 0..65535.
Pyy = 5000 * log10(Y.*conj(Y) * 10000);

%Шаг 9. Нарисуем график нашего спектра
%Подготовим ось X - ось частот от 0 до 8191 Гц
f =  (0:8191);
% и нарисуем график. Рисуем не весь массив, а только половину (вторая половина Pyy
является зеркальным отображением первой).
plot(f,Pyy(1:8192))
title ('Power spectral density')
xlabel('Frequency (Hz)')

%Шаг 10. Выведем наш спектр в файл на языке C, чтобы его можно было использовать
в микроконтроллерных системах.
stream_c = fopen('c:\spectrdata.c','wt');
fprintf (stream_c, '//таблица, сгенерированная MATLAB 6.5\n');
fprintf (stream_c, 'const unsigned short SpectrData [%i] = \n{\n', 8192);
for k=1:8192
   if (1==bitand(k, 7))
      fprintf(stream_c  , ' ');
   end;
   %пишем значение спектра в файл, отрицательные значения превращаем в 0
   spectr = round(Pyy(k));
   if (spectr>=0)
        fprintf(stream_c  , '%i', spectr);
   else
        fprintf(stream_c  , '%i', 0);
   end;
   if (not((k+1)==8193))
       fprintf(stream_c  , ',');
   end;
   if (0==bitand(k, 7))
       fprintf(stream_c  , '\n');
   end;
end
fprintf (stream_c, '};\n');
%файл spectrdata.c подготовлен и закрыт
fclose(stream_c);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

//таблица, сгенерированная MATLAB 6.5
const unsigned short SpectrData [8192] = 
{
    30792,30778,30785,30821,30848,30881,30976,30990,
    31129,31207,31304,31452,31548,31782,31931,32152,
    32397,32675,32968,33293,33608,33746,32872,34102,
    ...
    6132,8973,7742,12968,4600,12594,6098,10607,
    7250,12628,5761,11078,6584,10585,10397,2429,
    8683,12370,2227,8499,11136,5533,10788,11201
};