Программирование ARM MATLAB: как сгенерировать сигнал и его спектр Mon, November 20 2017  

Поделиться

нашли опечатку?

Пожалуйста, сообщите об этом - просто выделите ошибочное слово или фразу и нажмите Shift Enter.


MATLAB: как сгенерировать сигнал и его спектр Печать
Добавил(а) microsin   

На простом (относительно) примере показывается генерация сигнала и получение его спектра. Код взят из примера Help -> Demos -> MATLAB -> Mathematics -> FFT for Spectral Analysis и переработан под условия задачи.

Условие задачи: имеем ряд частот и соответствующих им амплитуд. Необходимо из этих частот сгенерировать сигнал в течение 1 секунды (оцифровка идет с частотой 16384 Гц) и построить его спектр в виде массива unsigned short SpectrData [8192] (каждый элемент массива 16-битный, всего 8192 значения). Номер элемента массива должен соответствовать частоте (от 0 до 8191 Гц), а значение элемента должно соответствовать логарифму амплитуды сигнала (логарифм нужен для лучшей наглядности представления малых амплитуд сигнала по сравнению с большими). Приведенный далее текст можно просто скопировать и вставить в MATLAB как код, сохранить как m-файл и запустить на выполнение.

 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);

В результате получим вот такую картинку спектра:

matlab-01-FFT-spectr.GIF

И такой текстовый файл spectrdata.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
};
 

Добавить комментарий


Защитный код
Обновить

Top of Page