Главная arrow Программирование arrow ARM arrow MATLAB: как сгенерировать сигнал и его спектр Saturday, June 24 2017  
ГлавнаяКонтактыАдминистрированиеПрограммированиеСсылки
UK-flag-ico.png English Version
GERMAN-flag-ico.png Die deutsche Version
map.gif карта сайта
нашли опечатку?

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

Поделиться:

MATLAB: как сгенерировать сигнал и его спектр Версия для печати
Написал microsin   
22.06.2009

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

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

Последнее обновление ( 22.06.2009 )
 

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

:D:lol::-);-)8):-|:-*:oops::sad::cry::o:-?:-x:eek::zzz:P:roll::sigh:

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

< Пред.   След. >

Top of Page
 
microsin © 2017