Преобразование Фурье и его свойства для ЦОС

Преобразование Фурье и его свойства для ЦОС 1) Создание приложения на формах Visual Studio 2017 (2019) C++ Поскольку все примеры курса далее приведены на Си (а не на C#) , приводится пример создания приложения на Forms на языке C++. Программы можно писать на C++ с оболочкой VS 3017 C++ или на C# , но тогда требуется переделывать синтаксис с Си на C#. Возможно также написание на Java, также с подключением бесплатных библиотек, например, графиков функций (в VS 2017 хватит chart). Ссылка на создание форм в VS 2017 C++ https://social.msdn.microsoft.com/Forums/vstudio/en-US/e6fbde42-d872-4ab3-8000-41ab22a4a584/visual-studio-2017-windows-forms?forum=winformsdesigner Основные шаги. 1. Установить, если не установлен C++/CLR 2. Создать пустой CLR проект 3. Проект - Свойства - Компоновщик - Система выбрать Windows (/SUBSYSTEM:WINDOWS) 4. Проект - Свойства - Компоновщик - Дополнительно - Точка входа набрать Main 5. Проект - Добавить новый элемент - UI - Форма Windows Form Эта ошибка будет повторяться, требуется открывать в конструкторе правой кнопкой в MyForm.h. 6. Скопировать в MyForm.cpp код: #include "MyForm.h" using namespace System; using namespace System::Windows::Forms; [STAThreadAttribute] void Main(array^ args) { Application::EnableVisualStyles(); Application::SetCompatibleTextRenderingDefault(false); Project1::MyForm form; Application::Run(%form); } При запуске: Далее как в C#. Вариант приложения Qt для C++ 2) Создание приложения на Qt языка C++ с использованием Visual Studio 2017 Установка Qt для работы с Visual Studio 2017 (2019) . https://medium.com/@paryleevatou/how-to-setting-qt-gui-using-c-with-visual-studio-2017-dbb50406f99 1. Установить Qt open source, подключить компоненты: MSVC2017, sources, charts, data visualization, debugger support 2. Установить Qt Visual Studio Tools https://marketplace.visualstudio.com/items?itemName=TheQtCompany.QtVisualStudioTools-19123 3. Добавить путь QT VS Tools - QT Options - Add Создание проекта Qt GUI через Qt Запуск F5 через Visual Studio 2017 Запуск F5 Визуальная работа с формами Открыть файл: QtGuiApplication1.ui Компоненты как в "чистом" Qt. 3) Сигналы, используемые в работе Для обработки с помощью преобразования Фурье представлены сигналы 4 типов, снятых у пациентов поликлиник медицинскими приборами в 90-х годах. Сигналы в файлах представлены в виде набора отсчетов либо по строкам, либо в виде одного столбца. 126 126 126 126 126 126 126 126 126 126 126 126 126 125 126 126 126 126 126 126 126 125 125 125 126 126 127 129 129 .......................... или 820 819 828 818 ...... Коды для кардио и рео сигналов находятся в диапазоне [0;255], вело и спиро сигналов диапазон [0;1023]. Примеры сигналов каждого типа, графики цифровых кодов сигналов. а) Cardio. Сигналы, полученные цифровым кардиографом. Пример графика отсчетов сигнала cardio28.txt Reo. Сигналы цифрового реографа. reo01.txt в кодах Velo. Сигналы кардиографа, пациент крутит педали велоэргометра. velo11.txt Spiro. Пациент дышит в трубку и измеряется объем выдыхаемого воздуха во времени. spiro9.txt Кроме медицинских сигналов в работе используются 3 стандартных тестовых сигнала: ступенька, пилообразный и треугольный сигналы. Значения отсчетов сигнала вещественные числа. Для всех трех типов сигналов задана его длина (во всех случаях N=720, объяснение далее). Для ступеньки Data[i]=1, i=N/2 Для пилообразного сигнала Data[i]= i * (1.0 / (N / 2.0)), i=N/2; Для треугольного сигнала Data[i] = i * (1.0 / (N / 2)), i=N/2 4) Основные понятия цифровых сигналов Основные задачи при переходе от аналогового сигнала к цифровому. Исходный аналоговый сигнал, например, с датчиков электрокардиографа, измеряющих напряжение, поступает на вход алфавитно-цифрового преобразователя (АЦП), который генерирует дискретные отсчеты сигнала через равные интервалы времени. Время между отсчетами называется интервалом выборки tв и измеряется в сек. Обратная величина называется частотой дискретизации сигнала fд=1/tв и измеряется в герцах (Гц). Например, если интервал выборки равен 1/100 сек, то частота дискретизации равна 100 Гц, т.е. в секунду генерируется 100 отсчетов. При этом частоту дискретизации стремятся сделать как можно больше, чтобы было как можно меньше потери информации. Пусть сигнал периодический с периодом T или с частотой fc=1/T. Тогда интервал выборки должен быть меньше полупериода сигнала или, иначе, частота дискретизации должна быть больше 2fc. В противном случае происходит потери информации о сигнале. Величина 2fc называется частотой Найквиста. Пример дискретизации синусоиды с частотой меньшей частоты Найквиста. У всех сигналов в данной лабораторной работе fд=360 Гц. 360 отсчетов соответствует одной секунде измерений. Для тестовых сигналов условно 720 отсчетов соответствует 2 сек, т.е. частота дискретизации 360 Гц. Другая проблема, это перевод аналогового значения сигнала в цифровое значение, которое называется уровнем квантования. Для этого настраивается вход АЦП на нижнее и верхнее значения сигнала. Нижнее значение соответствует выходу равному нулю, верхнее соответствует 2N, N – разрядность АЦП. На рис. N=3. Сигналы в вариантах для лабораторной работы получены либо на 8-разрядном АЦП(уровни 0-255), кардио и рео сигналы, либо 10-разрядном(уровни 0-1023) вело и спиро сигналы. 5) Единицы измерения цифровых сигналов в работе Кардио. Сигнал с датчика измеряется в милливольтах (мВ) положительное или отрицательное значение, диапазон уровней 0-255, 0 мВ соответствует коду 127, 60 уровней вверх (код 187) равно 1 мВ, соответственно 60 значений вниз равно -1 мВ. Рео. Сигнал измеряется в миллиомах (м Ом), емкостное сопротивление кожи. Диапазон 0-255, значения сопротивления только положительные. Нулевой уровень - нулевая проводимость. Каждые 100 значений добавляют добавляют 0,1 мОм. Вело. Сигнал также измеряется в мВ. Диапазон 0-1023. Нулевой уровень – 512 (0 мВ). Смещение на 240 меняет сигнал на 1 мВ. Спиро. Сигнал измеряется в литрах. Диапазон 0-1023. Нулевой уровень 512. +100 – вдох 1 литр, -100, соответственно, выдох 1 литр. В физических единицах измерения по оси OX откладывается время в сек. Для всех сигналов 360 отсчетов соответствует 1 сек. Для кардио сигналов отсчеты в мВ по оси OY пересчитываются: DataE[i]=(Data[i]-127)/60; i=0,N-1; пример графика cardio28.txt в физических единицах. Реосигнал DataE[i]=0.1*Data[i]/100; i=0,N-1; пример графика reo01.txt в физических единицах. Вело сигнал DataE[i]=(Data[i]-512)/240; i=0,N-1; пример графика velo11.txt в физических единицах. Спиро сигнал DataE[i]=(Data[i]-512)/100; i=0,N-1; пример графика spiro9.txt в физических единицах. В лабораторной работе следует использовать следующие параметры тестовых сигналов: частота дискретизации 360 Гц, длительность сигналов 2 сек (соответственно, количество отсчетов 720), нижняя граница 0, верхняя граница 1. Графики сигналов имеют вид: Ступенька Пила Треугольник 6) Разложение периодического сигнала в ряд Фурье Рассмотрим периодический непрерывный синусоидальный сигнал: Параметры сигнала. A – амплитуда сигнала измеряется в метрах (м) T – период сигнала (сек) w – угловая или круговая частота (рад./ сек) - фаза сигнала (радиан) Круговая частота w связана с обычной частотой fc, т.е. с числом повторений периода в одну секунду (измеряемая в Гц) соотношением: Например, если , то Гц, т.е. за 1 сек синусоида повторяется 2 раза. Соответственно период T=1/2 сек Основная теорема Фурье анализа. Любой непрерывный периодический сигнал s(t) с периодом T можно представить в виде бесконечной суммы: (1) где: - начальная частота Частота kw0 называется k- гармоникой представления ряда Фурье. Коэффициенты ak и bk рассчитываются по формулам: (2) a0 вычисляется по формуле для ak. - среднее значение сигнала за период T Пределы интегрирования можно сместить, часто используют период: [-T/2; T/2], что удобнее для интегрирования. Если s(t) четная функция, то bk=0, если нечетная, то ak=0 . В общем виде нужны все коэффициенты. Формулу (1) можно записать в виде только синусов или косинусов:( вещественная форма ряда Фурье), но тогда добавляются величины : (3) - эквивалентное представление периодического сигнала s(t). Связь между коэффициентами ak,bk формулы (1) и Ak, формулы (3) определяются выражениями: Представление сигнала в виде: позволяет придать смысл коэффициентам Ak, . Ak, k=0,1,2... – называется амплитудный спектр сигнала, он показывает величину вклада каждой гармоники. - фазовый спектр сигнала, на какой начальный угол смещена k-я гармоника. Для теоретического анализа ряд (1) или эквивалентный (3) можно записать также и в комплексной форме. Преобразование основано на формуле Эйлера для экспоненты с чисто мнимым показателем: где: j – мнимая единица. j2=-1, y - вещественное число. тогда cos (y) можно записать в виде: (4) Подставив (4) в вещественную форму ряда Фурье (1) получим: (5) представление сигнала s(t) в виде комплексных экспонент. Формула для расчета комплексных коэффициентов ck получается в виде комплексного интеграла: (6) Тогда получаются формулы связей (7): Ak=2|ck| амплитудный спектр фазовый спектр ak=2Re(ck) коэффициенты cos в разложении (1) bk=-2Im(ck) коэффициенты sin в разложении (1) Амплитудный спектр прямоугольного сигнала. Функция четная, в ряде Фурье (формула 1), остаются только косинусы. Для непрерывного сигнала ступеньки можно аналитически вычислить интегралы по формулам (2) , в результате получаем: Ak=ak= (8) График амплитудного спектра (первые 20 отсчетов Ak по формуле 8, всего для непрерывной ступеньки бесконечное количество отсчетов. рис 1 Амплитудный спектр непрерывного ступенчатого сигнала первые 20 отсчетов 7) Дискретное преобразование Фурье (ДПФ) Рассмотрим представление непрерывного периодического синала s(t) в виде ряда (см. 6 формула (5)) комплексные коэффициенты ck вычисляются по формуле: (6 формула (6)) T - период сигнала На практике имеется оцифрованный сигнал s(t) , т.е. N равноотстоящих отсчетов на интервале [0;T]. Обозначим эти отсчеты s0,s1,…,sN-1. При переходе к дискретным отсчетам сумма в бесконечных пределах формула (5) превращается в конечную сумму: m=0;N-1; (9) Комплексные коэффициенты ck для дискретных отсчетов рассчитываются заменой интегралов формулы (6) суммой (1/T переходит в 1/N): m=0;N-1; (10) Формулы (9) и (10) одинаковы, за исключением множителя 1/N. Формула (10) вычисления коэффициентов cm называется прямым дискретным преобразованием Фурье (ДПФ) сигнала s0,s1,…,sN-1. Формула (9) восстановления сигнала s0,s1,…,sN-1 по заданным комплексным коэффициентам c0,c1,…,cN-1 называется обратным дискретным преобразованием Фурье. Если отсчеты исходного сигнала sm - комплексные числа, формулы (9) (10) также справедливы. Комплексный дискретный сигнал рассматривается в основном в теоретических исследованиях. Если отсчеты sm вещественные, то вычисленные комплексные коэффициенты cm симметричны относительно N/2. cm=cN-m-1; m=0;N-1. На практике вычисляют раздельно вещественные и мнимые части коэффициентов ck разложения сигнала в ряд Фурье. ck=ak+jbk. Обозначим s[0],...,s[N-1] дискретные отсчеты вещественного сигнала, тогда вещественные и мнимые коэффициенты Фурье вычисляются по формулам, которые получаются из (9): , k=0,N-1; (11) , k=0,N-1 Обратная формула получения отсчетов сигнала sk по вещественным и мнимым частям коэффициентов разложения сигнала ak, bk имеет вид: (12) Т.к. формулы (11) и (12) практически одинаковы, вычисление прямого и обратного дискретного преобразования Фурье выполняется в одной функции с параметром t. Если t=0, вычисляется прямое преобразование, если t=1, обратное. Циклы вычисления ДПФ имеют вид. Программа 1 вычисления прямого и обратного ДПФ t=0 - прямое, t=1 обратное s[N] - массив отсчетов сигнала a[N],b[N] - массивы вещественной и мнимой части преобразования сигнала s[N] P=3.1415926; for(k=0;k f0 f1 001 - 100 -> f4 f2 010 - 010 -> f2 f3 011 - 110 -> f6 f4 100 - 001 -> f1 f5 101 - 101 -> f5 f6 110 - 011 -> f3 f7 111 - 111 -> f7 реверс бит в обратном порядке [f0,f1,...,fN-1] -> [f0,fi1,...,fN-1] Пусть f=[f0,f1,...,fN-1] массив, после перестановки элементов. Итерационный алгоритм (комплексные числа). Значения преобразования Фурье записываются в те же элементы массива f0,f1,...,fN-1 . N – длина вектора (степень двойки), P – показатель степени 2P=N. for(s=1;s<=P;s+) { m=2^s; wm=exp(2*Pi*j/m); // wm – комплексное j – мнимая единица for(k=0;k=0) { k=((I & Low)<>Shift); R=R | k; Shift -= 2; Low =Low <<1; High=High>>1; } return R; } тестовый пример: K=reverse(33,6); //43 ( 101011) , 2^6 -> 53 (110101) Код БПФ на языке Си: int i, M, k, j, N2, v, m; double Sa,fh; double WmR, WmI; double WR, WI; double tr, ti, ur, ui; double x, y; M = 1; N2 = N; while (N2 > 1) { N2 = N2 / 2; if (N2 > 1) M++; } if (t == 0) { for (i = 0; i < N; i++) { a[i] = s[i]; b[i] = 0; } } for (i = 0; i < N; i++) { k = reverse(i, M); if (i < k) { Sa = a[i]; a[i] = a[k]; a[k] = Sa; Sa = b[i]; b[i] = b[k]; b[k] = Sa; } } for (v = 1; v <= M; v++) { m = 1; for (j = 1; j <= v; j++) { m = m * 2; } WmR = cos(2 * 3.1415926 / m); WmI = sin(2 * 3.1415926 / m); for (k = 0; k < N; k += m) { WR = 1; WI = 0; for (j = 0; j <= m / 2 - 1; j++) { tr = a[k + j + m / 2] * WR - b[k + j + m / 2] * WI; ti = a[k + j + m / 2] * WI + b[k + j + m / 2] * WR; ur = a[k + j]; ui = b[k + j]; a[k + j] = tr + ur; b[k + j] = ti + ui; a[k + j + m / 2] = ur - tr; b[k + j + m / 2] = ui - ti; x = WR * WmR - WI * WmI; y = WR * WmI + WI * WmR; WR = x; WI = y; } } } Пример. spiro9.txt. Исходная длина сигнала N=6157. Чтобы использовать БПФ отбрасываем 2061 отсчет. N=4096. Амплитудный, фазовый спектры восстановленные сигналы полностью совпадают. Время прямого преобразования DPF = 6,45 сек Время прямого преобразования BPF = 0,004 сек. Быстрее в 1600 раз. 12) Ускоренное преобразование Фурье При работе с БПФ часто приходится либо отбрасывать существенную часть сигнала, либо добавлять отсчеты до степени двойки (копируя, например, фрагменты). Существуют ряд методов, которые медленнее БПФ, но быстрее ДПФ. Пусть N можно представить в виде произведения: N=L*M, при этом: L=2P . Чем больше L , тем ближе алгоритм к быстрому преобразованию Фурье и быстрее работает. Сначала применяется последовательно быстрое преобразование Фурье для L отчетов с шагом M. f0,f1,...,fN-1 – отчеты сигнала (или преобразованные Фурье, комплексные для выполнения обратного преобразования.) Вычисляется последовательно БПФ M раз для L отчетов сигнала с шагомм M: f0,fM, f2M...,f(L-1)M – БПФ с записью на те же места f1,fM+1, f2M+1...,f(L-1)M+1 – БПФ с записью на те же места ......................................... fM-1,fM-1+M, f2M+M-1...,f(L-1)M+M-1 – БПФ с записью на те же места Далее вычисляется комплексный массив F – преобразование Фурье(или раздельно действительные и мнимые части F1 и F2) for(s=0;s 1) { N2 = N2 / 2; if (N2 > 1) M++; } if (t == 0) { for (i = Z; i < N*H; i += H) { a[i] = Data[i]; b[i] = 0; } } for (i = 0; i < N; i++) { k = reverse(i, M); p = Z + i * H; q = Z + k * H; if (p < q) { Sa = a[p]; a[p] = a[q]; a[q] = Sa; Sa = b[p]; b[p] = b[q]; b[q] = Sa; } } for (s = 1; s <= M; s++) { m = 1; for (j = 1; j <= s; j++) { m = m * 2; } WmR = cos(2 * Pi / m); WmI = sin(2 * Pi / m); for (k = 0; k < N; k += m) { WR = 1; WI = 0; for (j = 0; j <= m / 2 - 1; j++) { p = k + j + m / 2; q = k + j; p = Z + p * H; q = Z + q * H; tr = a[p] * WR - b[p] * WI; ti = a[p] * WI + b[p] * WR; ur = a[q]; ui = b[q]; a[q] = tr + ur; b[q] = ti + ui; a[p] = ur - tr; b[p] = ui - ti; x = WR * WmR - WI * WmI; y = WR * WmI + WI * WmR; WR = x; WI = y; } } } } void BPFn(int t) { int k, i, P, L, M, N2 , r, m,v; double Pi; Pi = 3.1415026; N2 = N; L = 1; P = 0; k = 1; while (k == 1) { if ((N2 % 2) == 0) { N2 = N2 / 2; L = L * 2; P++; } else k = 0; } M = N2; for (i = 0; i < M; i++) { BPFh(0, i, M, L); // вызов } for (v = 0; v < M; v++) { for (r = 0; r < L; r++) { a1[r + v * L] = 0; b1[r + v * L] = 0; for (m = 0; m < M; m++) { a1[r + v * L] += a[m + r * M] * cos(2 * Pi*m*(r + v * L) / N) - b[m + r * M] * sin(2 * Pi*m*(r + v * L) / N); b1[r + v * L] += a[m + r * M] * sin(2 * Pi*m*(r + v * L) / N) + b[m + r * M] * cos(2 * Pi*m*(r + v * L) / N); } } } for (i = 0; i