В реальной жизни на динамические системы большое количество факторов оказывают влияние. Эти факторы чаще всего носят непостоянный характер, так что в результате некоторые особенности поведения системы непредсказуемы просто из-за невозможности учесть все взаимодействия. Но в теории колебаний и нелинейной динамике все же есть инструменты для исследования именно случайных процессов. Дело в том, что при невозможности предсказать кратковременные изменения системы мы зачастую можем оценить поведение усредненных статистических характеристик. Различные средние показатели, спектры, вероятностные распределения и другие характеристики могут сообщить много важного о модели и моделируемой системе в прошлом и в будущем, хотя не смогут помочь предсказать близжайшее изменение. Пишите в комментарии, если интересно подробнее узнать про стохастические модели, а здесь мы обсудим, как оценивать и анализировать источники равномерного белого некоррелированного шума. Рассмотренные здесь источники шума являются достаточно хорошими и поэтому характеристики будут близкими к теоретическим идеальным.
На рисунках 1 и 2 приведены реализации модельного равномерного шума из разных источников. Эти зависимости при наложении сильно друг от друга будут отличаться, но их характеристики практически одинаковы, так что ниже буду рассматривать только источник шума ran() из книги Numerical Recipes (рисунок 2). Алгоритм получения значений модельного шума взят из упомянутой книги и выглядит следующим образом (с некоторой моей адаптацией):
//Uniform noise generator////////////////////////
float ran (long *idum) {
static long MBIG = 1000000000;
static long MSEED = 161803398;
static long MZ = 0;
static float FAC = (1.0/MBIG);
static int inext,inextp;
static long ma[56];
static int iff=0;
long mj,mk;
int i,ii,k;
/*Инициализация*/
if (*idum < 0 || iff == 0) {
iff=1;
mj=labs(MSEED-labs(*idum));
mj %= MBIG;
ma[55]=mj;
mk=1;
for (i=1;i<=54;i++) {
ii=(21*i) % 55;
ma[ii]=mk;
mk=mj-mk;
if (mk < MZ) mk += MBIG;
mj=ma[ii];
}
for (k=1;k<=4;k++)
for (i=1;i<=55;i++) {
ma[i] -= ma[1+(i+30) % 55];
if (ma[i] < MZ) ma[i] += MBIG;
}
inext=0;
inextp=31;
*idum=1;
}
if (++inext == 56) inext=1;
if (++inextp == 56) inextp=1;
mj=ma[inext]-ma[inextp];
if (mj < MZ) mj += MBIG;
ma[inext]=mj;
return mj*FAC;
}
//Uniform noise generator////////////////////////
Для того, чтобы понять, является ли этот модельный источник на самом деле равномерным, нужно вычислить распределение вероятности (должно быть равномерным в некотором заданном отрезке). Спектр показывает, является ли шум белым (равномерное распределение на всех частотах) или цветным (распределение в спектре неравномерное). Поистине белым шум может быть только теоретически, а в реальности все источники шума ограничены по частоте сверху. Наконец, случайность проверяется корреляцией (автокорреляция должна начинаться с единицы и резко снижаться до нуля и больше не увеличиваться). Если корреляция вдруг увеличивается после нулевого смещения, шум является коррелированным и не вполне случайным.
На рисунке 3 приведено распределение для источника некоррелированного белого равномерного шума ran() в пределах от -10 до 10. Видно, что все значения случайной величины близки к одному и тому же значению 0.05. Кстати, если длину интервала (-10:10) умножить на это значение 0.05, получится единица, как и должно быть для распределений вероятности. Аналогичная зависимость наблюдается для источника rand() из базовых библиотек C/C++.
Спектр источника некоррелированного белого равномерного шума ran() (рисунок 4) также является, как и ожидается, равномерным. На графике нет выраженных максимумов и минимумов, как это было бы в случае источника цветного шума. Аналогичный спектр наблюдается для источника rand() из базовых библиотек C/C++.
Наконец, автокорреляция случайной величины источника некоррелированного белого равномерного шума ran() (рисунок 5) равна единице при отсутствии смещения (первая точка на графике) и близка к нулю во всех остальных точках. На графике показаны первая 1000 точек, а расчеты были сделаны для 10000 точек и после 1000 никаких новых трендов не обнаружено. Скорее всего, если продолжать считать корреляцию и дальше, один или второй источники шума могут оказаться скоррелированными с собой, но здесь мы говорим в первую очередь о том, как должны выглядеть характеристики. Для источника rand() из базовых библиотек C/C++ наблюдается аналогичная зависимость.
Таким образом, источники шума ran() и rand() являются источниками некоррелированного белого равномерного шума и вполне подходят для моделирования стохастических систем. В продолжение этой статьи построим аналогичные характеристики и проведем анализ для нормального шума и полетов Леви. Продолжение статьи о шуме.
p.s. Полную программу для исследования шума можете скачать на сайте кафедры радиофизики и нелинейной динамике СГУ.
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
Функция на C/C++ для вычисления характеристик шума
int uniform_distribution (long *idum, double dist_min, double dist_max, std::string tsname, std::string distname, std::string corrname, std::string specname) {
static int iter_num = 20000;
float *iterations; iterations = new float[iter_num]; int marker = 0;
static int corr_num = 10000;
float *meanx; meanx = new float[corr_num];
float *meany; meany = new float[corr_num];
float *mx; mx = new float[corr_num];
float *my; my = new float[corr_num];
float *cov; cov = new float[corr_num];
float *correlations; correlations = new float[corr_num];
int *corrcounters; corrcounters = new int[corr_num];
float epsilon = 1.0e-6; int counter_prec = 0; int prec_limit = 10;
float min, max; static int dist_num = 100;
int * distribution; distribution = new int[dist_num];
int dist_counter = 0; double dp = 0.0; double dmin = 0.0;
int spectra_num = 5000; int spectra_count = 0;
float *spectra_part_im; spectra_part_im = new float[spectra_num];
float *spectra_part_re; spectra_part_re = new float[spectra_num];
float *spectra_im; spectra_im = new float[spectra_num];
float *spectra_re; spectra_re = new float[spectra_num];
bool generate = true;
int counter = 0;
while (generate) {
iterations[marker] = ran(idum)*(dist_max-dist_min)+dist_min;
if (counter==(iter_num-1)) {
for (int i=0;i<corr_num;i++)
corrcounters[i]=0;
for (int i=0;i<corr_num;i++) {
meany[i] = meanx[i] = 0.0;
for (int j=i;j<iter_num;j++) {
meanx[i] += iterations[j-i];
meany[i] += iterations[j];
corrcounters[i]++;
}
meanx[i]/=(1.0*corrcounters[i]);
meany[i]/=(1.0*corrcounters[i]);
mx[i] = my[i] = cov[i] = 0.0;
for (int j=i;j<iter_num;j++) {
float diffx = (iterations[j-i]-meanx[i]);
float diffy = (iterations[j]-meany[i]);
mx[i] += diffx*diffx;
my[i] += diffy*diffy;
cov[i] += diffx*diffy;
}
}
min = iterations[0];
max = iterations[0];
for (int j=1;j<iter_num;j++) {
if (iterations[j]<min)
min = iterations[j];
if (iterations[j]>max)
max = iterations[j];
}
for (int j=0;j<dist_num;j++) {
distribution[j] = 0;
}
dmin = min;
double dmaxloc = max;
if ((dmin<0.0)&&(dmaxloc>0.0)) {
if (fabs(dmin)>fabs(dmaxloc)) dmaxloc = -dmin;
else dmin = -dmaxloc;
}
dp = (dmaxloc-dmin)/dist_num;
for (int j=0;j<iter_num;j++) {
dist_counter++;
for (int i=0;i<dist_num;i++) {
if ((iterations[j]>(dmin+dp*i)) && (iterations[j]<(dmin+dp*(i+1))))
distribution[i]++;
}
}
for (int k=0; k<spectra_num; k++) {
spectra_re[k] = 0.0;
spectra_im[k] = 0.0;
for (int j=0;j<iter_num;j++) {
spectra_re[k] += iterations[j] * cos(2.0 * M_PI * k * j / iter_num);
spectra_im[k] -= iterations[j] * sin(2.0 * M_PI * k * j / iter_num);
}
spectra_re[k] /= iter_num;
spectra_im[k] /= iter_num;
}
spectra_count++;
}
else {
if (marker==(iter_num-1)) {
for (int k=0; k<spectra_num; k++) {
spectra_part_re[k] = 0.0;
spectra_part_im[k] = 0.0;
for (int j=0;j<iter_num;j++) {
spectra_part_re[k] += iterations[j] * cos(2.0 * M_PI * k * j / iter_num);
spectra_part_im[k] -= iterations[j] * sin(2.0 * M_PI * k * j / iter_num);
}
spectra_re[k]+=spectra_part_re[k]/iter_num;
spectra_im[k]+=spectra_part_im[k]/iter_num;
}
spectra_count++;
}
}
if (counter>=iter_num) {
bool prec_now = true;
for (int i=0;i<corr_num;i++) {
int delay = marker - i;
if (delay<0) delay+=iter_num;
float oldmeanx = meanx[i];
float olddiffx = (iterations[delay]-oldmeanx);
float oldmeany = meany[i];
float olddiffy = (iterations[marker]-oldmeany);
meanx[i]=oldmeanx+olddiffx/(corrcounters[i]+1.0);
meany[i]=oldmeany+olddiffy/(corrcounters[i]+1.0);
mx[i] += olddiffx*(iterations[delay]-meanx[i]);
my[i] += olddiffy*(iterations[marker]-meany[i]);
cov[i] += olddiffx*olddiffy*(corrcounters[i]+0.0)/(corrcounters[i]+1.0);
float corr_new = cov[i]/(sqrt(mx[i]*my[i]));
if (fabs(corr_new-correlations[i])>epsilon) prec_now = false;
correlations[i] = corr_new;
corrcounters[i]++;
}
if (prec_now) {
counter_prec++;
if (counter_prec==prec_limit)
generate = false;
}
else
counter_prec = 0;
if (iterations[marker]<min)
min = iterations[marker];
if (iterations[marker]>max)
max = iterations[marker];
dist_counter++;
for (int i=0;i<dist_num;i++) {
if ((iterations[marker]>(dmin+dp*i)) && (iterations[marker]<(dmin+dp*(i+1))))
distribution[i]++;
}
}
if (marker==0) std::cout << "finished iteration numuber: " << counter << "\n";
counter++; marker++; if (marker==iter_num) marker = 0;
}
std::ofstream ts(tsname);
ts << "# min=" << min << " max=" << max << "\n";
for (int j=0;j<iter_num;j++) {
int cm = marker+j; if (cm>=iter_num) cm-=iter_num;
ts << counter-iter_num+j+1 << " " << iterations[cm] << "\n";
}
ts.close();
std::ofstream dist(distname);
dist << "# min=" << min << " max=" << max << "\n";
for (int i=0;i<dist_num;i++) {
dist << i << " " << dmin+dp*i << " " << dmin+dp*(i+1) << " " << 1.0*distribution[i]/(dist_counter*dp) << " " << distribution[i] << "\n";
}
dist.close();
float mincorr = correlations[1];
float maxcorr = correlations[1];
for (int i=2;i<corr_num;i++) {
if (correlations[i]<mincorr)
mincorr=correlations[i];
if (correlations[i]>maxcorr)
maxcorr=correlations[i];
}
std::ofstream corr(corrname);
corr << "# mincorr=" << mincorr << " maxcorr=" << maxcorr << "\n";
for (int i=0;i<corr_num;i++) {
corr << i << " " << meanx[i] << " " << meany[i] << " " << mx[i]/corrcounters[i] << " " << my[i]/corrcounters[i] << " " << cov[i]/corrcounters[i] << " " << correlations[i] << "\n";
}
corr.close();
std::ofstream spec(specname);
for (int k=0; k<spectra_num; k++) {
spectra_re[k]/=spectra_count;
spectra_im[k]/=spectra_count;
spec << 1.0*k/iter_num << " " << spectra_re[k] << " " << spectra_im[k] << " " << sqrt(spectra_re[k]*spectra_re[k]+spectra_im[k]*spectra_im[k]) << " " << atan(spectra_im[k]/spectra_re[k]) << "\n";
}
spec.close();
delete [] spectra_re;
delete [] spectra_im;
delete [] spectra_part_re;
delete [] spectra_part_im;
delete [] distribution;
delete [] correlations;
delete [] cov;
delete [] my;
delete [] mx;
delete [] meany;
delete [] meanx;
delete [] iterations;
return 0;
}
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
