Scilab - генератор сигналов.

Недавно возникла необходимость в программном многофункциональном звуковом генераторе, позволяющим менять форму сигнала, частоту, воспроизводить различные типы псевдошумовых сигналов. Конечно, существует достаточное количество коммерческих и свободных программных продуктов, но как-то не интересно было копаться в них, тем более, если есть возможность самому сделать такую программку, и достаточно быстро, используя возможности Scilab. Несколько слов, что такое Scilab? Scilab (www.scilab.org) - это некоммерческая система компьютерной математики, предназначенная для численных вычислений. В этом пакете реализованы алгоритмы численного решения для многих типов задач: решение нелинейные уравнений и систем, решение задач линейной алгебры, решение задач оптимизации, дифференцирование и интегрирование, задачи обработка экспериментальных данных, решение обыкновенных дифференциальных уравнений и систем и т.д.. Scilab позволяет создавать и редактировать различные типы графиков от1d до 3d. Кроме этого, в состав Scilab входит Хcos — система визуального компьютерного моделирования, аналогичная Simulink (Matlab). Возможности Scilab могут быть расширены внешними программами и модулями, написанными на разных языках программирования. Для ОС Linux существует также возможность интеграции в Scilab системы символьных вычислений Maxima (http://maxima.sourceforge.net/ru/).

Scilab , как и программы Matlab, Оctave, позволяет работать со звуковой картой компьютера, записывать и воспроизводить звуковые файлы, чем и воспользуемся. Обсудим, как можно генерировать различные типы сигналов. Стоит отметить, что в Scilab очень просто оперировать с любыми элементарными (в том числе тригонометрическими) функциями. Будем использовать функцию sin(x) для генерации звука нужной нам частоты. Аргумент х означает какое-то абстрактное значение угла, пока не установлена частота выборки (сэмплинг) значений функции в единицу времени (секунду) и частота звука. Учитывая это, запишем значение аргумента х, используя синтаксис Scilab:

x= 2*%pi*f*t ;

где- %pi – значение пи, f- частота звука [Гц], t-время звучания [сек]. Как известно, современные звуковые карты компьютеров поддерживают частоты дискретизации 8,0, 11,025, 22,05 и 44,1 кГц, а профессиональные и больше. Частота дискретизации (или частота семплирования, англ. sample rate) это - частота взятия отсчетов непрерывного во времени сигнала при его дискретизации. Измеряется в герцах. Частота дискретизации (оцифровки) сигнала должна быть как минимум в два раза больше максимальной частоты входного сигнала (теорема Котельникова-Найквиста). Если человеческая речь, например, занимает полосу частот до 3-4 кГц, то для ее оцифровки потребуется частота 8 кГц и т.п. Таким образом, нужно связать длительность звучания t с частотой дискретизации звуковой карты Fs. Можно записать: t=N/Fs; где N-общее количество отсчётов аргумента (пропорционально длительности нашего сигнала). Тогда аргумент х будет выглядеть следующим образом:

x= 2*%pi*(f/Fs)*N;

Теперь получим для примера синусоидальный звуковой сигнал длительностью 2 секунды частотой 440 Гц. Пусть частота дискретизации звуковой карты 44100 Гц. Запишем

N=0:88200;// длительность сигнала в отсчётах

Fs=44100;// частота дискретизации

f=440;// частота сигнала

x =2*%pi *(f/Fs)*N;// аргумент функции

y=sin(x);//синусоидальный сигнал

sound(y)// воспроизведение сигнала

plot(x,y)//график сигнала

plot(x(1:500),y(1:500))//график пятьсот первых точек сигнала

                                                Рис1. Графики синусоидального сигнала

В Scilab для воспроизведения звука есть две команды: playsnd и sound (см. help).
Если мы зададим наш сигнал в виде y = А*sin(x), где А –число, то можем легко изменять амплитуду сигнала. Стоит отметить также, что в Scilab есть команда

soundsec(t[,Fs])

с помощью которой можно сгенерировать вектор сэмплов для сигнала длительностью t сек. при частоте дискретизации Fs [Гц]. Т.е. если t = soundsec(2,Fs), то можно записать наш аргумент:

    х=2*%pi*f*t.

В случае необходимости гармоник, их можно задать следующим образом:

          y2 = sin(2*x), y3 = sin(3*x), y4 = sin(4*x)…

Давайте посмотрим, как можно задать сигналы другой формы и записать их в файл с расширением .wav. Пусть частота для всех форм сигналов будет 440 Гц.

Для синусоидального сигнала это будет выглядеть так:

f=440; // частота сигнала
sin_s = sin(2*%pi*f*t);
xx=pwd();//текущая директория
sinfile = xx+'\sin'+string(f)+'.wav';
wavwrite(sin_s, Fs, bits, sinfile);
где bits может быть: 8, 16, 24, 32 бита – количество битов, используемых для записи информации для каждого сэмпла (16 – по умолчанию).

 

 

Пилообразный сигнал:

pila=2*(f*t*floor(0.5+f*t));
xx=pwd();//текущая директория
pila_file = xx+'\pila'+string(f)+'.wav';
wavwrite(pila, Fs, bits, pila_file);

Рис.2 Пилообразный сигнал

 

 

Треугольный
triangle=(2/%pi)*asin(sin(2*%pi*f*t));
xx=pwd();//текущая директория
tri_file = xx+'\triangle'+string(f)+'.wav';
wavwrite(triangle, Fs, bits, tri_file);

 

Рис.3 Треугольный сигнал

 

Прямоугольный
square=sign(sin(2*%pi*f*t));
xx=pwd();//текущая директория
sq_file = xx+'\square'+string(f)+'.wav';
wavwrite(square, Fs, bits, sq_file);

Рис.4 Прямоугольный сигнал

 

Полученные звуковые файлы можно прочитать в Scilab, используя команду [y,Fs,bits]=wavread(имя файла) (подробности см. в справке).  
А теперь посмотрим, как можно сгенерировать белый гауссов  шум и добавить его в сигнал. Белый шум — стационарный шум, спектральные составляющие которого равномерно распределены по всему диапазону задействованных частот, с дельта-функцией Дирака в качестве автокорреляционной функции. Примеры белого шума: белый свет (откуда пошло название шума), шум близкого водопада, дробовой шум на клеммах большого сопротивления, шум стабилитрона (при протекании очень малого тока).   В Scilab, белый гауссов шум можно легко сгенерировать, используя генератор случайных чисел rand или grand. При этом не будем думать об автокорреляционной функции.  Используем для этой цели команду rand(N,'normal'), которая сгенерирует случайные числа с нормальной плотностью распределения со средним m=0, и среднеквадратичным отклонением sigma=1. Добавим эту команду в наш код:

 

N=0:88200;// длительность сигнала в отсчётах
Fs=44100;// частота дискретизации
f=440;// частота сигнала
x =2*%pi *(f/Fs)*N;// аргумент функции
y=sin(x);//синусоидальный сигнал
y1=rand(N,'normal');//гауссов шум
sound(y)// воспроизведение сигнала
sound(y+y1);
plot(x,y)//график сигнала
z=y+y1;//добавление шума к сигналу
plot(x(1:500),z(1:500))//график- пятьсот первых точек сигнала

Рис.5 Синусоидальный сигнал с добавленным шумом, белый дискретный и белый гауссов шум в переменных: амплитуда, время

 

На рис.5 можно увидеть, как выглядят белый и белый гауссов шум в переменных: амплитуда, время и синусоидальный сигнал с добавленным шумом.

И напоследок, рассмотрим, как сгенерировать "розовый" (pink), мерцательный (фликкер) шум. Например, так шумит отдаленный водопад, так как высокочастотные составляющие звука затухают в воздухе сильнее низкочастотных. Спектральная плотность мощности "розового" шума определяется как 1/f (плотность обратно пропорциональна частоте) и равномерно убывает в логарифмической шкале частот. Спектральная плотность такого сигнала по сравнению с белым шумом затухает на 3 децибела на каждую октаву. Послушать, как звучат цветные шумы можно по этой ссылке: https://ru.wikipedia.org/wiki/%D0%91%D0%B5%D0%BB%D1%8B%D0%B9_%D1%88%D1%83%D0%BC
            Посмотрим, как можно сгенерировать розовый шум, используя возможности Scilab.  Отметим, что розовый шум можно получить  фильтрацией белого шума или исходя из определения. Рассмотрим оба случая. В первом случае, был адаптирован под Scilab код (для Matlab и Octave), : https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html

Код для Scilab генератора розового шума:

 

Nx = 2^16;  //number of samples to synthesize
B = [0.049922035 -0.095993537 0.050612699 -0.004408786];
A = [1 -2.494956002   2.017265875  -0.522189400];
AA=flipdim(A,2)
z=poly(AA,'y','coeff');
nT60 = round(log(1000)/(1-max(abs(roots(z))))); //T60 est.
v = rand(1,Nx+nT60,'normal');// Gaussian white noise: N(0,1)
x = filter(B,A,v);    // Apply 1/F roll-off to PSD
end=length(x);
x = x(nT60+1:end);   // Skip transient response
sound(x)

Исходя из определения свойств "розового" шума, можно адаптировать под  Scilab следующий код (https://www.mathworks.com/matlabcentral/fileexchange/42919-pink--red--blue-and-violet-noise-generation-with-matlab-implementation/content/pinknoise.m):

 

// y - row vector of pink (flicker) noise samples
N=2^16;//N - number of samples
rem=N-fix(N/2)*2;
if rem==1
M = N+1;
else
M = N;
end
// generate white noise
x = rand(1:M,'normal');
// FFT
X = fft(x);// fast Fourier transform (FFT)
// prepare a vector for 1/f multiplication
NumUniquePts = M/2 + 1;
n = 1:NumUniquePts;
n = sqrt(n);
// multiplicate the left half of the spectrum so the power spectral density
// is proportional to the frequency by factor 1/f, i.e. the
// amplitudes are proportional to 1/sqrt(f)
X(1:NumUniquePts) = X(1:NumUniquePts)./n;

// prepare a right half of the spectrum - a copy of the left one,
// except the DC component and Nyquist frequency - they are unique
X(NumUniquePts+1:M) = real(X(M/2:-1:2)) -%i*imag(X(M/2:-1:2));

// IFFT
y = ifft(X);//inverse fast Fourier transform
// prepare output vector y
y = real(y(1, 1:N));
// ensure unity standard deviation and zero mean value
y = y - mean(y);
yrms = sqrt(mean(y.^2));
y = y/yrms;
sound(y)

Рис.6 Розовый шум

 

Сравнение качества генераторов шума это отдельный вопрос, к которому мы вернемся в дальнейшем.