Спектральная Плотность Мощности - $5 за пару строк кода

[vi19540315]

В приложенном архивном файле проект PSD на C#.
Проект читает файл rr.txt и пишет файл pp.txt
Входной файл это временной ряд из 300 чисел зафиксированных в течение 339095 миллисекунд.
Каждое число ряда это время в миллисекундах прошедшее после предыдущего измерения.
Программа считывает временной ряд в массив и применяет к нему дискретное преобразование Фурье.
Результаты преобразования выводятся в файл pp.txt

В программе не хватает нескольких строк для расчета мощности каждой гармоники
и вместо значения мощности выводятся нули.

Нужно добавить код для расчета мощности и вывода ее в составе строки файла.

Расчитанная для каждой гармоники мощность не должна отличаться от приведеной в файле ee.txt
более чем на 2%.

Для расчета мощности можно написать либо оригинальный код либо найти готовую функцию
в какой-либо библиотеке.

Стоимость $5.
Оплата через PayPal после проверки работы кода на моем компьютере.

[WorldMaster]

Фигня какая то. А где данные сигнала то??

Сообщение от vi19540315

Входной файл это временной ряд из 300 чисел зафиксированных в течение 339095 миллисекунд.

Тут написано что файл содержит сигнал с периодом 1130,31 мс.

Сообщение от vi19540315

Каждое число ряда это время в миллисекундах прошедшее после предыдущего измерения.

А тут уже получается что в файле время в мс.


Так где измерения то?? где сам сигнал?

Более того. В конечном файле у вас единицы должны быть в ms в квадрате на герц а мощность получается в герцах.
Для перевода еще необходимо понимать физику процесса. нужен коэффициент, который зависит от физических свойств сигнала.

[jillitil]

Там всё не то.

Что такое БПФ я знаю.

♦Именно поэтому функции скормили временные отрезки вместо амплитудных.
♦Замеры через разные промежутки времени, чего быть не должно.
♦На выходе получили комплексное число, а надо досчитать до скаляра. ТС скопипастил код не разобравшись.
♦Желание получить на выходе физическую величину, при этом нет входных данных.
♦Причём в "эталонном" ее.тхт есть намёк что расчёт возможен.

В программе не хватает нескольких строк

Сумма всех значений в спектрограмме и есть искомая величина.
Чтоб до этого добраться надо писать весь код.

не должна отличаться более чем на 2%.

Без оконной функции (как в коде) результат = рюмка плюс-минус ведро.


ТС, не подскажите что значит квадратная миллисекунда на герц? Стандартно PSD ед.измерения это Вт/Гц или дБм/Гц.


В целом скопипащенная ТС-ом и доработанная функция кажется работает:

Type Complex
    R As Double
    I As Double
End Type

Dim F(MaxN) As Complex, FFT(MaxN) As Double

Sub DFT()
    Maxi = 0
    For I = 1 To MaxN \ 2
        F(I).R = 0
        F(I).I = 0
        For J = 1 To MaxN
                R1 = N(J).R
                I1 = N(J).I
                R2 = Cos(w * (J - 1) * (I - 1))
                I2 = -Sin(w * (J - 1) * (I - 1))
                R3 = R1 * R2 - I1 * I2
                I3 = R1 * I2 + I1 * R2                
                F(I).R = F(I).R + R3
                F(I).I = F(I).I + I3
        Next
        FFT(I) = Sqr(F(I).R ^ 2 + F(I).I ^ 2)
        If FFT(I) > Maxi Then Maxi = FFT(I)
    Next
    
' normalize + PSD
    For I = 1 To MaxN \ 2
        FFT(I) = FFT(I) / Maxi
        PSD = PSD + FFT(I)
    Next I
End Sub

[vi19540315]

В описании задачи указано что проект читает файл rr.txt
Файл rr.txt приложен к теме
В нем перечень чисел
Именно они и составляют ряд исходных данных
Временные замеры между интервалами действительно не одинаковы но этим можно пренебречь
Считаем что одинаковые (ну для применения ДПФ)
Данные ряда это интервалы между ударами сердца
Так оно бьется - не совсем равномерно
Мощность в миллисукундах в квадрате на герц
Шумы уже устранены и грубые ошибки импульсов в сигнале тоже
Детрендинг не нужен
Применение оконной функции приветствуется
Я про нее смутное представление имею
Если я не на все ответил или еще нужна информация - пишите

В конце нужно примерное совпадение результата расчета мощности с указанной в файле ee.txt
(А там точно оконная функция применялась)

[vi19540315]

Похоже тут для оценки мощности применялся метод периодограмм Велча и оконная функция Хэмминга. Так что парой строк кода тут явно не обойдешься.
Поэтому тема закрывается, а господин jillitil если хочет то может сообщить мне свой PayPal акк для получения $5.

Но все таки странно. Если есть Фурье разложение определенного временного ряда то мощность спектра одной определенной гармоники заданной частоты должна бы иметь одно и то же значение хоть ты тресни. А тут применяют какие-то периодограммы и функции разные и получают разные результаты и что - все они правильные? Похоже как я купил что-то за 10 копеек а мне с сотни дают любую сумму сдачи и всякий вариант правильный? Не понимаю такой математики.

[jillitil]

Сообщение от vi19540315

Но все таки странно.... должна бы иметь одно и то же значение хоть ты тресни.

Набор данных во времени это сигнал. Сигнал можно охарактеризовать, например макс-мин, RMS, пик-пик, среднее арифметическое, и т.п.

Самый первый и примитивный графический вид сигнала это осциллограмма или амплитудо-временная характеристика (АВХ) [моя отсебятина]. Иное представление этого же сигнала амплитудо-частотная хар-ка (АЧХ)

БПФ переводит АВХ в АЧХ. Но тут есть проблема: фаза и частота в функции на практике не совпадают с таковыми в обрабатываемом сигнале. В результате преобразования на АЧХ появятся несуществующие частоты (вы их некорректно назвали гармониками). Называется Spectral leakage (утечка спектра) см. Рис и сравните с предыдущим две высшие частоты. Функция дискретно проходит, а в данных я забил дробную частоту, они размазались. Для уменьшения их появления применяется "окно".

Окно Хемминга(Hamming): хорошо убирает несуществующие, но может подавить присутствующую частоту.
Окно Ганна (Hann): теоретически хуже Хемминга, но не давит существующие частоты полностью. Лучше пользоваться им в решениях, где важно ничего не упустить, ценой лишнего шума.

Сообщение от vi19540315

все они правильные?

Абсолютно правильных решений как в школе нет. Задачи сильно сложные.

[vi19540315]

Спасибо за разъяснение. Картина существенно прояснилась.
Предложение насчет пейпала остается в силе.

[vi19540315]

Остается еще один не понятный для меня вопрос.

Метод Велча включает следующие шаги:
- Разделите исходный сигнал на перекрывающиеся или неперекрывающиеся сегменты равной длины.
- Примените оконную функцию к каждому сегменту для уменьшения эффекта утечки спектра.
- Вычислите периодограмму для каждого сегмента, которая определяется как квадрат модуля дискретного преобразования Фурье (DFT) сегмента.
- Усредните периодограммы по всем сегментам для получения оценки PSD.

Значит имеем входной ряд из 300 элементов и продолжительностью сьемки 340 секунд.
Допустим я разбил свой входной ряд из 300 элементов на 59 участков по 10 элементов от 1 по 10, от 6 по 15, от 21 по 30 и т.д.
Применил к каждому участку оконную функцию Хемминга.
Применил к каждому участку ДПФ.
На каждом участке получил по 10 магнитуд.

-----

Нашел формулу усреднения периодограмм в общем виде

$$ \text{PSD}(f) = \frac{1}{N}\sum_{i=1}^{N}\frac{1}{M }\left|\text{FFT}(x_i)\right|^2 $$

где:

$\text{PSD}(f)$ - оценка спектральной плотности мощности (PSD) на частоте $f$;
$N$ - число сегментов;
$M$ - длина каждого сегмента;
$x_i$ - $i$-й сегмент сигнала;
$\text{FFT}(x_i)$ - преобразование Фурье (FFT) $i$-го сегмента;
$|\text{FFT}(x_i)|^2$ - квадрат амплитуды спектра $i$-го сегмента.
Таким образом, мы берем среднее значение квадрата амплитуды спектра для каждого сегмента и делим его на общее число сегментов,
чтобы получить усредненную оценку PSD на данной частоте $f$.

Заметим, что PSD - это функция частоты, поэтому для оценки PSD на всех возможных частотах,
необходимо применять эту формулу для каждой из частот, на которых мы хотим оценить PSD.

-----

Вопрос:
Можно описать эту формулу нормальным человеческим языком чоб можно было это записать в виде программного кода?
Что с чем складывать и на сколько делить чтобы получить 300 значений Спектральной Оценки Мощности для исходного ряда из 300 элементов?

------

Убрав служебные разметочные символы получаем

PSD(f) = (1 / N) * Summa(po i ot 1 do N) ( (1 / M) * |FFT(x_i)|^2 )

где:

PSD(f) - оценка спектральной плотности мощности (PSD) на частоте f;
N - число сегментов;
M - длина каждого сегмента;
x_i - i-й сегмент сигнала;
FFT(x_i) - преобразование Фурье (FFT) i-го сегмента;
|FFT(x_i)|^2 - квадрат амплитуды спектра i-го сегмента.
Таким образом, мы берем среднее значение квадрата амплитуды спектра для каждого сегмента и делим его на общее число сегментов,
чтобы получить усредненную оценку PSD на данной частоте f.

Получается что частота f присутствует в правой части через x_i ?
И я не понимаю - откуда наберется 300 штук PSD(f)

Так чтоб сопоставить их с ДПФ исодного ряда.

-----

Частота PSD связана с номером элемента в сегменте через формулу:

f_i = i * (F_s / N)

где f_i - частота i-го элемента в сегменте,
F_s - частота дискретизации сигнала,
а N - количество элементов в сегменте.

Номер сегмента не имеет непосредственной связи с частотой PSD.
Вместо этого, номер сегмента используется для усреднения спектров.
Обычно сегменты перекрываются на определенном проценте, например, 50%.
Для каждого сегмента вычисляется спектр, затем спектры усредняются.

Т.е. на выходе получаем 10 значений PSD т.к. в сегменте было 10 элементов, а сколько было сегментов не имеет значения.
Но мне же нужно 300 значений PSD. И где они?