Берем сигнал длительностью
Если есть возможность выбрать
Далее, создаем массив комплексных чисел, длиной
ShortComplex arr[N];
Теперь заполняем его. Допустим исходный сигнал у нас в массиве
for(int i= 0; i < N; ++i)
{
arr[i].re= x[i];
arr[i].im= 0.0;
}
Теперь осталось применить алгоритм БПФ. Если
fft(arr, k, false);
Если же мы не можем гарантировать, что
universal_fft(arr, N, false);
После исполнения этой функции массив
Для обратного преобразования надо вызвать ту же функцию с последним аргументом
universal_fft(arr, N, true);
for(int i= 0; i < N; ++i)
x[i]= arr[i].re;
Полученный результат представляет собой последовательность комплексных чисел в форме пар: (реальная часть, мнимая часть). Но для понимания физической структуры сигнала желательно преобразовать их в амплитуды, частоты и фазы гармоник.
Для преобразования можно использовать формулы:
Самый первый элемент с индексом
Часто бывает дано не суммарное время сигнала
В сумме постоянная составляющая и все гармоники должны дать исходный сигнал.
Также по ходу дела можно исправить зеркальный эффект. Если
взять сигнал, который состоит из единственной гармоники с частотой
Для устранения зеркального эффекта можно просто обнулить вторую половину спектра с индексами после
В процессе преобразования могут возникнуть ошибки округления, из-за которых появятся ненулевые элементы там, где должен быть ноль. Если мы хотим устранить этот эффект, нужно задать некоторый нижний предел и отбрасывать гармоники с амплитудами меньше этого предела.
И, наконец, последний штрих: функция синуса которая начинается от нуля, может оказаться
нагляднее, чем функция косинуса. Поэтому есть смысл представить гармоники в виде
синусов:
В результате алгоритм получения наглядного представления будет таким:
//Это комплексные числа, результат прямого дискретного преобразования Фурье
//ShortComplex - это структура с полями: double re, im;
ShortComplex arr[N];
//это индекс комплексного числа в массиве arr
int i;
//это частота дискретизации
double nSamplesPerSec;
//убираем зеркальный эффект, просто отбрасывая вторую половину
int Nmax= (N + 1) / 2;
//мы хотим получить массив гармоник.
//это массив амплитуд, массив частот и массив фаз для каждой гармоники
double *freq= new double[Nmax];
double *amp= new double[Nmax];
double *phase= new double[Nmax];
//это индекс гармоники. в конце алгоритма он будет равен количеству найденных гармоник
int j= 0;
//это нижний предел амплитуды гармоники. тут может быть число больше нуля
//в зависимости от того, какие амплитуды мы считаем ничтожно малыми.
double limit= 0.001;
//вспомогательная переменная для оптимизации
double abs2min= limit * limit * N * N;
//получаем постоянную составляющую
if (arr[i]->re >= limit)
{
amp[j]= arr[i].re / N;
freq[j]= 0.0;
phase[j]= 0.0;
++j;
}
++i;
//получаем остальные гармоники
for(i= 1; i < Nmax; ++i)
{
double re= arr[i].re;
double im= arr[i].im;
//это квадрат модуля комплексного числа arr[i]
double abs2= re * re + im * im;
//отбрасываем слишком слабые гармоники
if (abs2 < abs2min)
continue;
//вычисляем апмлитуду. 2.0 - для устранения зеркального эффекта
amp[j]= 2.0 * sqrt(abs2) / N;
//вычисляем фазу косинуса в радианах
phase[j]= atan2(im, re);
//преобразуем косинус в синус. M_PI2 = пи/2, M_PI = пи
//в результате фаза будет в диапазоне от -пи/2 до +пи/2
phase[j]+= M_PI2;
if (phase[j] > M_PI)
phase[j]-= M_2PI;
//можно еще преобразовать радианы в градусы
phase[j]= phase[j] * 180.0 / M_PI;
//получаем частоту
freq[j]= (nSamplesPerSec * i) / N;
++j;
}
Получив список частот, фаз и амплитуд, можно воспользоваться какой-нибудь программой, чтобы построить график или таблицу. Например, в очередной версии программы Bard, которую я сейчас пишу, это выглядит так:
Если выполнить прямое, а потом обратное преобразование Фурье, то должны получиться те же самые числа с небольшими отклонениями из-за неизбежных ошибок округления. В частности, если исходная последовательность была действительной (мнимые части равны нулю), то после прямого преобразования Фурье могут появиться ненулевые мнимые части, но после обратного преобразования они снова обнулятся.
Однако, если после прямого преобразования внести изменения в массив и потом выполнить обратное преобразование, то могут появиться элементы с ненулевыми мнимыми частями. Если вы применяете преобразования Фурье для обработки действительных сигналов, эти ненулевые мнимые части надо просто игнорировать.
В частности, если вы исправляете зеркальный эффект, вы (после прямого преобразование) обнуляете элементы
с индексами выше
Данный эффект подробно расматривался ранее, а здесь приведен среди других эффектов.
Если частота исходной гармоники не равна в точности дной из частот
Здесь зеленым показана амплитуда спектральных гармоник, фиолетовым - фаза (как видите,
в районе нужной частоты она делает резкий скачок). Исходной гармоникой здесь было
колебание с амплитудой 30000, фазой 0 и частотой 440,2 Гц. Время дискретизации
В англоязычной литературе этот эффект называется smearing (размазывание) или leakage (рассеяние). Это связано с тем, что мощность исходной гармоники как бы рассеивается. Если вы посмотрите на предыдущую иллюстрацию, то увидите, что высота пика не дотягивает до правильной отметки в 30000, но вместе с дополнительными линиями, окружающими пик, мощность получается такой, как нужно.
А вот эта картинка -
- получена совсем другим способом. Здесь была гармоника 440 Гц,
но ее амплитуда в течение секунды линейно затухала от 30000 до 0. Картинки очень похожи,
так что пойди разбери - то ли причина размазывания в том, что гармоника не попала точно
в величину
Различие там есть, просто оно не видно на картинке. При непопадании на
На левой картинке - размазывание от несовпадения частот, на правой - от затухания. Так что в целом для различения этих случаев можно применять критерий: если фаза максимальной гармоники где-то посредине между фазами соседних гармоник, то это - модуляция, а если фаза максимальной гармоники гораздо близка к одной из "соседок", то это - несовпадение частот.
Вот еще картинка:
- здесь амплитуда линейно затухала от 30000 до 15000.
- здесь амплитуда нелинейно нарастала от 0 до 30000.
- эта скромная картинка получена с добавлением пульсации: амплитуда синусоидально увеличивается-уменьшается 10 раз за секунду. Две дополнительные частоты по бокам от основной обеспечили пульсацию основной частоты.
- а здесь все то же самое, но пульсация не 10 раз в секунду, а чуть-чуть медленнее. В секунду не уложилось целое число пульсаций, так что возник эффект размазывания вокруг дополнительных частот.
- наконец, вот этот хаос возникает в результате "вибрато" - то есть, периодического изменения не амплитуды, а частоты. В данном случае на периоде 0.8 секунд произошло 8 периодических изменений частоты на четверть тона.
Релятивисты и позитивисты утверждают, что "мысленный эксперимент" весьма полезный интрумент для проверки теорий (также возникающих в нашем уме) на непротиворечивость. В этом они обманывают людей, так как любая проверка может осуществляться только независимым от объекта проверки источником. Сам заявитель гипотезы не может быть проверкой своего же заявления, так как причина самого этого заявления есть отсутствие видимых для заявителя противоречий в заявлении.
Это мы видим на примере СТО и ОТО, превратившихся в своеобразный вид религии, управляющей наукой и общественным мнением. Никакое количество фактов, противоречащих им, не может преодолеть формулу Эйнштейна: "Если факт не соответствует теории - измените факт" (В другом варианте " - Факт не соответствует теории? - Тем хуже для факта").
Максимально, на что может претендовать "мысленный эксперимент" - это только на внутреннюю непротиворечивость гипотезы в рамках собственной, часто отнюдь не истинной логики заявителя. Соответсвие практике это не проверяет. Настоящая проверка может состояться только в действительном физическом эксперименте.
Эксперимент на то и эксперимент, что он есть не изощрение мысли, а проверка мысли. Непротиворечивая внутри себя мысль не может сама себя проверить. Это доказано Куртом Гёделем.
Понятие "мысленный эксперимент" придумано специально спекулянтами - релятивистами для шулерской подмены реальной проверки мысли на практике (эксперимента) своим "честным словом". Подробнее читайте в FAQ по эфирной физике.
|
|
 |
