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

 В первой и второй части (part1, part2) рассматривалось, как с помощью   Scilab и  Хcos  можно создать генераторы сигналов и шумов. Но хотелось бы понять, насколько эти генераторы правильно воспроизводят заданные сигналы. Поэтому, в этой части займёмся спектральным анализом полученных сигналов. В Scilab встроено достаточно большое количество функций, позволяющих проводить обработку сигналов https://www.scilab.org/scilab/features/scilab/signal_processing.
 Воспользуемся этими возможностями. Рассмотрим, вначале, насколько правильно генерируется синусоидальный звуковой сигнал и сравним его с сигналом, записанным в файл sin440.wav(см.part1). Также, к коду для генерации синусоидального сигнала добавим код для фурье-преобразования и затем сравним ширину пиков на полувыоте:

N=0:88200;// длительность сигнала в отсчётах
Fs=44100;// частота дискретизации
f=440;// частота сигнала
x =2*%pi *(f/Fs)*N;// аргумент функции
y=sin(x);//синусоидальный сигнал
Y=fft(y,-1);
Ymag=abs(Y);
N=max(N);
Ymag=Ymag(1:N/2);
F = 0:((N/2)-1);
F = F * Fs/N;
subplot (2,1,1)
plot2d(F(860:900),Ymag(860:900));
plot(F(860:900),22050*ones(1,41),'--k')
xtitle('|X(f)|','f (Hz)', ' ')
Сигнал из файла:
[y,Fs,bits]=wavread('D:\scilab\sin440.wav');Fs,bits
[a,N]=size(y);t=0:(N-1);
t=t/Fs;
Y=fft(y,-1);
Ymag=abs(Y);
N=max(N);
Ymag=Ymag(1:N/2);
F = 0:((N/2)-1);
F = F * Fs/N;
subplot (2,1,2)
plot(F(860:900),Ymag(860:900));
plot(F(860:900),22050*ones(1,41),'--k')

Полученные спектры представлены на рисунке 1. Как видно, ширина пика и положение для обоих случаев одинаковое.

Рис.1

Теперь рассмотрим частотные спектры наших генераторов шума.  В качестве реперных сигналов шума возьмём файлы с белым и розовым шумами с сайта: https://ru.wikipedia.org/wiki/Цвета_шума  
Длительность этих сигналов 10 сек. Обрежем их до 2 секунд. Кстати, в Scilab легко можно обрезать длительность звуковых файлов (единственно, возникает проблема нехватки памяти, если файл длинный, но в версии Scilab 6.0 проблему устранили).

[y,Fs,bits]=wavread('D:\scilab\White_noise.wav');Fs,bits
N=2*Fs;// длительность сигнала 2 сек.
t=0:(N-1);
t=t/Fs;
y1=rand(1:N,'uniform')//генерация белого шума в Scilab
Y1=fft(y1,-1);//фурье-преобразование
Y=fft(y,-1);
Ymag=abs(Y);
Y1mag=abs(Y1);
N=max(N);
Ymag=Ymag(1:N/2);
Y1mag=Y1mag(1:N/2);
F = 0:((N/2)-1);
F = F * Fs/N;
subplot(2,1,1)
plot2d(F,Ymag);
xtitle('|X(f)|','f (Hz)', ' ')
subplot(2,1,2)
plot2d(F,Y1mag);
xtitle('|X(f)|','f (Hz)', ' ')
variance(y);variance(y1)//дисперсия сигнала из файла и из Scilab

Дисперсия сигнала из файла  0.008, для программного генератора rand 0.08. Как видно, дисперсии различаются на порядок. Если взять генератор белого гауссова шума rand( ,'normal'), то дисперсии окажутся примерно одинаковы. Поэтому, есть вопрос, каким образом формировался сигнал в файле. На рисунке 2 представлены спектры белого шума из файла и сгенерированного в Scilab.

Рис.2

Посмотрим теперь на частотные спектры "розового" шума, полученные в Scilab и реперном файле. Воспользуемся вышеприведённым кодом для "белого шума" и прочитаем файл для "розового" шума. В итоге, получим частотный спектр "розового" шума. На рисунке 3 –прямоугольники синего цвета. В случае программного генератора, код которого представлен ниже, спектр – звёздочки красного цвета. Как видно, поведение обоих спектров одинаковое.

N=2^16;//N - количество сэмплов
rem=N-fix(N/2)*2;
if rem==1
M = N+1;
else
M = N;
end
// генерация белого шума
x = rand(1:M,'normal');
// FFT
X = fft(x);// быстрое Фурье преобразование(FFT)
// подготовка вектора для умножения на 1/f
NUP = M/2 + 1;
n = 1:NUP;
n = sqrt(n);
X(1:NUP) = X(1:NUP)./n;
X(NUP+1:M) = real(X(M/2:-1:2)) -%i*imag(X(M/2:-1:2));
// IFFT
y = ifft(X);//обратное быстрое Фурье преобразование
// подготовка выходного вектора y
y = real(y(1, 1:N));
// обеспечивание равенства единицы стандартного отклонения и 0 =среднего значения
y = y - mean(y);
yrms = sqrt(mean(y.^2));
y = y/yrms;
Y=fft(y,-1);
Y1mag=abs(Y);
N=max(N);
Ymag=Y1mag(1:N/2);
F = 0:((N/2)-1);
F = F * Fs/N;
Ymag=10*log10(Ymag);
plot(F,Ymag,'*r');
xtitle('|X(f)|','f (Hz)', ' ')

Рис.3

В part1 рассматривались два алгоритма генерации розового шума. Давайте сравним их спектры. На рисунке 4 представлены спектры для обоих генераторов. Видно, что поведение спектров одинаковое, но один спектр смещён относительного другого. Синим цветом отображается спектр генератора, который рассматривали выше. Красным – спектр генератора в алгоритме, которого используется фильтрация белого шума.    

Рис.4

И, напоследок, сравним спектры "розового" шума, полученные с помощью алгоритма, исходящего из определения розового шума, и с помощью Xcos (генератор шума на рис.5)

 

Рис.5

На рисунке 6 видно, что спектры для двух генераторов практически совпадают.

Рис.6