Меня тут немного покритиковали за мою предыдущую статью. Я сделал генератор псевдослучайных чисел в ПЛИС платы Марсоход3бис и сказал, что мол визуально шумит - значит работает. Визуально распределение равномерное - значит - работает.
Но нет. Мне Alexandr говорит, что нужно найти периодичность, спектр, автокорреляцию и только потом что-то утверждать. Чтобы хоть как-то оправдаться, решил хотя бы найти и показать спектр шумового сигнала, полученного из моего генератора LFSR.
Кто не знает, что такое спектр, пожалуй дам ссылку на свою же очень давнюю статью.
Теперь к делу..
Программа для отрисовки спектра не совсем моя и не совсем написана мною. Скажем так: мы высококвалифицированные программисты, поэтому в конце концов нагуглили решение и создали свое из разных источников.
Для использования этой программы на питоне потребуется не только сам питон, но и многие дополнительные его библиотеки, например, такая как scipy. SciPy - это основанная на питоне экосистема для математики, научных и инженерных расчетов.
SciPy (pronounced “Sigh Pie”) is a Python-based ecosystem of open-source software for mathematics, science, and engineering.
Чтобы не мучаться с поиском и установкой всяких питоновских библиотек типа SciPy рекомендую устанавливать не чистый питон, а питон вместе со всеми этими пакетами: я использую WinPython.
Программа для рисования исходного сигнала и его амплитудного спектра вот:
import numpy as np from pylab import * from scipy import * import os def read_bin_file(filename): f=open(filename, 'rb') values = np.fromfile(f,dtype="int8") return values def plotSpectrum(y,Fs): n = len(y) # length of the signal k = arange(n) T = n/Fs frq = k/T # two sides frequency range frq = frq[range(n/2)] # one side frequency range Y = fft(y)/n # fft computing and normalization Y = Y[range(n/2)] plot(frq,abs(Y),'r') # plotting the spectrum xlabel('Freq (Hz)') ylabel('|Y(freq)|') Fs=1024.0 y = read_bin_file("samples_imp.bin") #y = read_bin_file("samples_saw.bin") #y = read_bin_file("samples_sin.bin") #y = read_bin_file("samples_rnd.bin") t = arange(0,len(y),1) # time vector subplot(2,1,1) plot(t,y) xlabel('Time') ylabel('Amplitude') subplot(2,1,2) plotSpectrum(y,Fs) show()
Эта программа читает в массив элементов бинарный файл, в котором каждый байт (как знаковое число) представляет из себя одну выборку сигнала. Какой именно файл читать задается просто в коде вызовом функции read_bin_file(). В программе несколько строк с read_bin_file(), но только одна работает - остальные закомментированы. Выбирайте тот файл, что вам интересен. Потом программа рисует сам сигнал во времени и его частотный спектр.
Чтобы убедиться, что программа рисования спектра правильно работает и все правильно рисует я дополнительно написал программу на C/C++, которая создает три двоичных файла с выборками для сигнала синусоиды, пилы и прямоугольных импульсов. Все три файла длиной 1024 байта. Вот программа, создающая тестовые файлы:
//
// test_waves.cpp : Defines the entry point for the console application. // #include "stdafx.h" #define _USE_MATH_DEFINES #include <math.h> #define NUM_SAMPLES 1024 int _tmain(int argc, _TCHAR* argv[]) { char arr_sin[NUM_SAMPLES]; char arr_saw[NUM_SAMPLES]; char arr_imp[NUM_SAMPLES]; double sin_step=2*M_PI*32/NUM_SAMPLES; for(int i=0; i<NUM_SAMPLES; i++) { arr_saw[i]=(i&0x3F)*4-127; arr_imp[i]=(i&0x40) ? 127 : -127; arr_sin[i]=127*sin(sin_step*i); } FILE * pFile; errno_t err = fopen_s(&pFile,"samples_sin.bin","wb"); fwrite(arr_sin,1,NUM_SAMPLES,pFile); fclose (pFile); err = fopen_s(&pFile,"samples_saw.bin","wb"); fwrite(arr_saw,1,NUM_SAMPLES,pFile); fclose (pFile); err = fopen_s(&pFile,"samples_imp.bin","wb"); fwrite(arr_imp,1,NUM_SAMPLES,pFile); fclose (pFile); return 0; }
Для тестовых сигналов пилы, меандра и синусоиды я умышленно сделал периоды сигналов кратными этому магическому числу 1024, так чтобы четное число периодов тестового сигнала четко укладывались по всему интервалу из 1024 выборок. Это сделано для того, чтобы спектр тестовых сигналов получился "чистым и красивым". Дело в том, что при вычислении спектра предполагается, что весь интервал выборок будет повторен бесконечное число раз влево и вправо. При склейке интервалов не должно быть швов, разрывов в сигнале. Тогда спектральные линии будут узкими и ярко выраженными. В программе на питоне я задаю частоту выборок в 1024 Гц. То есть весь интервал времени из тестовых файлов длиной 1024 байта - это одна секунда.
Итак, после запуска сишной программы у меня получились три тестовых файла samples_sin.bin, samples_saw.bin, samples_imp.bin. Посмотрим, как их отрисовывает питоновская программа:
Синусоида:
Видна одна четкая спектральная линия на частоте 32Гц. Так и должна быть, ведь в тестовом сигнале 1024 выборки на одну секунду как раз 32 периода синусоиды.
Посмотрим на прямогольные импульсы, меандр:
В тестовом периоде за одну секунду 8 импульсов и поэтому первая гармоника 8Гц. Потом видны, как и положено, нечетные третья гармоника 24Гц и пятая гармоника 40Гц.. Ну и так далее.
Пила:
Шестнадцать импульсов пилы за тестовый период в одну секунду. Первая гармоника 16Гц и потом идут четные гармоники 32Гц, 48Гц и так далее.. Похоже на правду.
Теперь, когда я уверен, что программа рисует спектр правильно можно попробовать посмотреть спектр сигнала из моего генератора псевдослучайной последовательности.
Я немного модифицировал программу приема случайных чисел из платы Марсоход3бис. Теперь программа не только принимает и показывает распределение чисел, но и имеет возможность записи выборок в файл.
В программе появилась кнопка - если ее нажать, то 1024 байт из псевдослучайной последовательности последовательности будут записаны в файл samples_rnd.bin.
Полученный файл можно передать на отрисовку питоновской программе.
Что получается? Вот что:
Не знаю.. кому как, а я бы сказал, что спектр сигнала довольно шумоподобный и довольно равномерный.
Амплидуда низкочастотных составляющих как-то странно ниже - возможно это признак того, что тестовый интервал в мои 1024 выборки плохо склеивается? Вот честно затрудняюсь сказать..
Ну вот такой получился результат.
Все программы в виде исходников, а так же тестовые сигналы в виде двоичных файлов можно взять здесь:
Подробнее...