Программа дискретного преобразования фурье

Программа дискретного преобразования фурье

Даётся программный код для прямого и обратного преобразования Фурье. Рассматривается быстрое преобразование Фурье.

Для расчёта преобразования Фурье обычно используется ускоренная процедура расчёта – т.н. быстрое преобразование Фурье (БПФ ). Это позволяет в значительной мере сократить процессорное время на достаточно сложные и ресурсоёмкие математические расчёты.

1 Комплексныечисла

Для начала нам потребуется вспомогательный класс, который будет описывать комплексные числа. Комплексные числа – это особый вид чисел в математике. Каждое комплексное число состоит из двух частей – действительной и мнимой. Сейчас нам достаточно знать о комплексных числах применительно к ДПФ то, что действительная часть комплексного числа хранит информацию об амплитуде сигнала, а мнимая – о фазе.

Код класса для описания комплексных чисел (разворачивается)

2 Прямое дискретное быстрое преобразование Фурье

Теперь приведу код функции, которая реализует расчёт прямого дискретного преобразования Фурье. Здесь встречается термин «бабочка». Не пугайтесь, это общепринятое название одной из элементарных операций, которые входят в алгоритм БПФ. Подробнее с термином «бабочка» можно ознакомиться здесь.

На вход функции передаётся массив комплексных чисел. Действительная часть которого представляет произвольный дискретный сигнал, с отсчётами через равные промежутки времени. Мнимая часть содержит нули. Число отсчётов в сигнале должно равняться степени двойки. Если ваш сигнал короче, то дополните его нулями до числа, кратного степени 2: 256, 512, 1024 и т.д. Чем длиннее сигнал, тем у рассчитанного спектра будет выше разрешение по частоте.

Код для расчёта прямого быстрого преобразования Фурье на VB.NET (разворачивается)

3 Обратное дискретное быстрое преобразование Фурье

Код для расчёта обратного быстрого преобразования Фурье на VB.NET (разворачивается)

Ну и конечно же, опишем использовавшийся метод, который проверяет число элементов переданного массива:

4 Проверка прямого и обратного преобразования Фурье

Теперь давайте проверим, что наши функции работают. Для этого пропустим произвольный сигнал через механизм прямого преобразования Фурье, а затем «соберём» его обратно с помощью обратного преобразования Фурье. Восстановленный сигнал должен практически совпадать с исходным. Ошибки округления, возникающие при работе с числами в компьютере, имеют место быть, поэтому сигналы не будут идентичны полностью, но их отклонение друг от друга должно быть пренебрежимо малым.

Для примера в качестве исходного сигнала возьмём функцию синуса и сформируем данные длиной 128 отсчётов вот таким образом:

Получим вот такой сигнал:

Исходный сигнал во временной области до преобразования Фурье

Здесь по оси X – номера отсчётов во временной области, по оси Y – амплитуда. Обратим внимание, что сигнал состоит только из действительных частей, а мнимая часть на всём отрезке равна "0".

Теперь передадим этот сигнал на вход функции FFT(). По полученным в ходе прямого преобразования Фурье массивам комплексных чисел построим два графика – действительной (Re) и мнимой (Im) частей спектра:

Действительная и мнимая части сигнала в частотной области

Здесь по оси X – отсчёты в частотной области, по оси Y – амплитуда. Чтобы получить реальные значения частоты, необходимо рассчитать их, учитывая, что "0" оси Y соответствует нулевой частоте, максимум оси Y соответствует частоте дискретизации.

Полученный спектр сигнала передадим функции обратного преобразования Фурье IFFT(). Получим массив комплексных чисел, где действительная часть будет содержать восстановленный сигнал:

Действительная и мнимая части восстановленного с помощью обратного преобразования Фурье сигнала

Как видно, восстановленный сигнал полностью повторяет исходный.

Преобразование Фурье

Определение преобразования Фурье

Преобразование Фурье является представлением изображения как сумма комплексных экспоненциалов переменных значений, частот и фаз. Преобразование Фурье проигрывает решающую роль в широком диапазоне приложений для обработки изображений, включая улучшение, анализ, восстановление и сжатие.

Если f ( m, n) является функцией двух дискретных пространственных переменных m и n , то двумерное преобразование Фурье f ( m, n) задано отношением

F ( ω 1 , ω 2 ) = ∑ m = − ∞ ∞ ∑ n = − ∞ ∞ f ( m , n ) e − j ω 1 m e − j ω 2 n .

Читайте также:  Как определить наименьший положительный период

Переменные ω 1 и ω 2 являются переменными частоты; их модули являются радианами на выборку. F ( ω 1, ω 2) часто называется представлением частотного диапазона f ( m, n) . F ( ω 1, ω 2) является комплексной функцией, которая является периодической и в ω 1 и в ω 2 с периодом 2 π . Из-за периодичности, обычно только область значений − π ≤ ω 1 , ω 2 ≤ π отображен. Обратите внимание на то, что F (0,0) является суммой всех значений f ( m, n) . Поэтому F (0,0) часто называется постоянным компонентом или компонентом DC преобразования Фурье. (DC обозначает постоянный ток; это — электротехническое понятие, которое относится к источнику питания постоянного напряжения, в противоположность источнику питания, напряжение которого отличается синусоидально.)

Инверсия преобразования является операцией, которая, когда выполняется на преобразованном изображении производит оригинальное изображение. Обратным двумерным преобразованием Фурье дают

f ( m , n ) = 1 4 π 2 ∫ ω 1 = − π π ∫ ω 2 = − π π F ( ω 1 , ω 2 ) e j ω 1 m e j ω 2 n d ω 1 d ω 2 .

Примерно разговор, это уравнение означает, что f ( m, n) может быть представлен как сумма бесконечного числа комплексных экспоненциалов (синусоиды) с различными частотами. Значение и фаза вклада на частотах ( ω 1, ω 2) даны F ( ω 1, ω 2) .

Визуализация преобразования Фурье

Чтобы проиллюстрировать, рассмотрите функциональный f ( m, n) , который еще равняется 1 в прямоугольной области и 0 везде. Чтобы упростить схему, f ( m, n) показывается непрерывной функцией, даже при том, что переменные m и n дискретны.

Прямоугольная функция

Следующие данные показывают, как сетчатый график, значение преобразования Фурье,

из прямоугольной функции, показанной в предыдущей фигуре. Сетчатый график значения является распространенным способом визуализировать преобразование Фурье.

Изображение значения прямоугольной функции

Пиком в центре графика является F (0,0) , который является суммой всех значений в f ( m, n) . График также показывает, что F ( ω 1, ω 2) имеет больше энергии на высоких горизонтальных частотах, чем в высоких частотах кадровой развертки. Это отражает то, что горизонтальные сечения f ( m, n) являются узкими импульсами, в то время как вертикальные сечения являются кадровыми импульсами. Узкие импульсы имеют больше высокочастотного содержимого, чем кадровые импульсы.

Другой распространенный способ визуализировать преобразование Фурье состоит в том, чтобы отобразиться

журнал | F ( ω 1 , ω 2 ) |

как изображение, как показано.

Журнал преобразования Фурье прямоугольной функции

Используя логарифм помогает произвести детали преобразования Фурье в областях, где F ( ω 1, ω 2) очень близко к 0.

Примеры преобразования Фурье для других простых форм показывают ниже.

Преобразования Фурье некоторых простых форм

Дискретное преобразование Фурье

Работа с преобразованием Фурье на компьютере обычно включает форму преобразования, известного как дискретное преобразование Фурье (DFT). Дискретное преобразование является преобразованием, значения ввода и вывода которого являются дискретными выборками, делая его удобным для манипуляции с компьютерами. Существует две основных причины использования этой формы преобразования:

Ввод и вывод ДПФ оба дискретен, который делает его удобным для манипуляций с компьютерами.

Существует алгоритм FAST для вычисления ДПФ, известного как быстрое преобразование Фурье (FFT).

ДПФ обычно задается для дискретного функционального f ( m, n) , который является ненулевым только по конечной области 0 ≤ m ≤ M − 1 и 0 ≤ n ≤ N − 1 . Двумерным ДПФ M на n и обратными отношениями ДПФ M на n дают

F ( p , q ) = ∑ m = 0 M − 1 ∑ n = 0 N − 1 f ( m , n ) e − j 2 π p m / M e − j 2 π q n / N p = 0 , 1 , . , M − 1 q = 0 , 1 , . , N − 1

f ( m , n ) = 1 M N ∑ p = 0 M − 1 ∑ q = 0 N − 1 F ( p , q ) e j 2 π p m / M e j 2 π q n / N m = 0 , 1 , . , M − 1 n = 0 , 1 , . , N − 1

Значения F ( p, q) коэффициенты ДПФ f ( m, n) . Коэффициент нулевой частоты, F (0,0) , часто называется "компонентом DC". DC является электротехническим понятием, которое обозначает постоянный ток. (Обратите внимание на то, что матричные индексы в MATLAB ® всегда запускаются в 1, а не 0; поэтому, элементы матрицы f (1,1) и F (1,1) соответствуют математическим количествам f (0,0) и F (0,0) , соответственно.)

Функции MATLAB fft , fft2 и fftn реализуют алгоритм быстрого преобразования Фурье для вычисления одномерного ДПФ, двумерного ДПФ и N-мерного ДПФ, соответственно. Функции ifft , ifft2 и ifftn вычисляют обратного ДПФ.

Отношение к преобразованию Фурье

Коэффициенты ДПФ F ( p, q) выборки преобразования Фурье F ( ω 1, ω 2) .

Читайте также:  Как пополнить баланс с мтс на билайн

F ( p , q ) = F ( ω 1 , ω 2 ) | ω 1 = 2 π p / M ω 2 = 2 π q / N p = 0 , 1 , . , M − 1 q = 0 , 1 , . , N − 1

Визуализация дискретного преобразования Фурье

Создайте матричный f , который подобен функциональному f ( m, n) в примере в Определении преобразования Фурье. Помните, что f ( m, n) равен 1 в прямоугольной области и 0 в другом месте. Используйте двухуровневое изображение, чтобы представлять f ( m, n) .

Вычислите и визуализируйте 30 30 ДПФ f с этими командами.

Дискретное преобразование Фурье, вычисленное без дополнения

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

Чтобы получить более прекрасную выборку преобразования Фурье, добавьте нулевое дополнение в f при вычислении его ДПФ. Нулевое дополнение и вычисление ДПФ могут быть выполнены на одном шаге с этой командой.

Эта команда нулевые клавиатуры f , чтобы 256 256 прежде вычислить ДПФ.

Дискретное преобразование Фурье, вычисленное с дополнением

Коэффициент нулевой частоты, однако, все еще отображен в левом верхнем углу, а не центре. Можно решить эту проблему при помощи функционального fftshift , который подкачивает квадранты F так, чтобы коэффициент нулевой частоты был в центре.

Получившийся график идентичен один показанный в Визуализации преобразования Фурье.

Приложения преобразования Фурье

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

Частотная характеристика линейных фильтров

Преобразование Фурье импульсного ответа линейного фильтра дает частотную характеристику фильтра. Функциональный freqz2 вычисляет и отображает частотную характеристику фильтра. Частотная характеристика Гауссова ядра свертки показывает, что этот фильтр передает низкие частоты и ослабляет высокие частоты.

Частотная характеристика гауссова фильтра

Смотрите Линейные фильтры Проекта в Частотном диапазоне для получения дополнительной информации о линейной фильтрации, отфильтруйте проект и частотные характеристики.

Выполните быструю свертку Используя преобразование Фурье

Этот пример показывает, как выполнить быструю свертку двух матриц с помощью преобразования Фурье. Ключевое свойство преобразования Фурье состоит в том, что умножение двух преобразований Фурье соответствует свертке связанных пространственных функций. Это свойство, вместе с быстрым преобразованием Фурье, формирует основание для быстрого алгоритма свертки.

Примечание: основанный на БПФ метод свертки чаще всего используется для больших входных параметров. Для маленьких входных параметров это обычно быстрее, чтобы использовать функцию imfilter .

Создайте две простых матрицы, A и B . A является матрицей M на n, и B является P-by-Q матрицей.

A нулевой клавиатуры и B так, чтобы они были, по крайней мере (M+P-1) (N+Q-1). (Часто A и B дополнены нулем к размеру, который является степенью 2, потому что fft2 является самым быстрым для этих размеров.) Пример заполняет матрицы, чтобы быть 8 8.

Вычислите двумерного ДПФ A и B с помощью функции fft2 . Умножьте эти двух ДПФ вместе и вычислите обратного двумерного ДПФ результата с помощью функции ifft2 .

Извлеките ненулевой фрагмент результата и удалите мнимую часть, вызванную ошибкой округления.

В статье «Моделирование РЗА. Часть первая.» нам удалось записать файл в формате Comtrade, содержащий токи режима КЗ, воспроизведенного на тестовой модели энергосистемы в среде MATLAB Simulink.

Далее перейдем к созданию цифровой модели устройства РЗА, которая позволит нам протестировать поведение защиты в данном режиме КЗ и других режимах работы энергосистемы.

Ограничим задачу моделирования рамками максимальной токовой защиты, упрощенный алгоритм функционирования которой приведен на рисунке 1.

Рисунок 1. Алгоритм максимальной токовой защиты

В рамках этой работы нам потребуется решить три основные задачи:

— выполнить цифровую обработку сигналов осциллограммы для получения интегральных параметров – действующих значений токов основной частоты 50 Гц;

— разработать и реализовать алгоритм пускового органа максимальной токовой защиты, который выполняет сравнение интегральных параметров сигналов с уставкой срабатывания;

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

Читайте также:  3055 Сканирование windows 7

Дискретная форма представления сигналов

Сигналы токов, полученные с модели MATLAB Simulink, представлены в дискретном виде с постоянной частотой дискретизации – последовательность мгновенных значений сигналов, зафиксированных через равные интервалы времени. Подобное представление имеют сигналы токов и напряжений на выходе аналого-цифровых преобразователей (АЦП) устройств РЗА.

Пример дискретного сигнала приведен на рисунке 2. Пунктирной кривой зеленого цвета показан исходный непрерывный сигнал, окружностями красного цвета – его дискретный эквивалент – мгновенные значения («A», «A1» и «A2»), зафиксированные на интервале времени от «t» (текущий момент времени) до «t2» с постоянным шагом дискретизации «TS». Величина, обратная шагу дискретизации «TS», носит название частоты дискретизации. Обозначим ее «FS».

Рисунок 2. Пример дискретного сигнала

В файле Comtrade, полученном с тестовой модели, шаг дискретизации «TS» составляет 500 мкс, что соответствует частоте дискретизации «FS», равной 2000 Гц, и количеству выборок на период основной частоты сигнала «NS», равному 40.

Цифровая обработка сигналов

Работа алгоритма максимальной токовой защиты, как правило, основана на сравнении действующего значения первой гармонической составляющей измеренного тока с уставкой срабатывания.

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

Текст m-функции в программе MATLAB, реализующей вычисление действующего значения с помощью дискретного преобразования Фурье приведен ниже:

function [Ak] = fourier_calc(A, k)

global Ns;

global n;

ak=2/Ns*sum(A.*cos(k*2*pi.*n/Ns));

bk=-2/Ns*sum(A.*sin(k*2*pi.*n/Ns));

Ak=((ak^2+bk^2)/2)^0.5;

end

Расчет действующего значения первой гармоники с помощью данной функции будем производить на каждом шаге дискретизации при получении каждой новой выборки, используя вновь полученное значение выборки и 39 значений выборок сигнала (Ns-1), полученных на предыдущих шагах дискретизации.

Ниже представлен текст «Protection.m» файла-сценария в программе MATLAB, который обеспечивает:

  • запуск на расчет модели энергосистемы в MATLAB Simulink и выгрузку результатов расчета режима КЗ в рабочую область программы;
  • вычисление действующих значений токов;
  • запись в формат Comtrade мгновенных значений токов, полученных из MATLAB Simulink, и действующих значений, рассчитанных в данном файле.

clear all

close all

%%%%%%%%% Исходные данные %%%%%%%%%%

global Ns, Ns=40; %количество выборок на период сигнала — глобальная переменная (используется в м-функциях и Simulink-моделях)

global n, n=0:Ns-1; %порядковый номер выборки в анализируемом окне — глобальная переменная (используется в м-функциях)

Моделирование КЗ

sim(‘ M.slx ‘) %запуск Simulink-модели на расчет

for m=1:fix(length(DATA)) %Вычисления выполняем на каждом шаге получения нового значения выборок сигналов токов

if m>=Ns %имеющееся количество выборок достаточно для выполнения расчета ДПФ

Ia=(DATA(m:-1:m-(Ns-1),1))’;

Ib=(DATA(m:-1:m-(Ns-1),2))’;

Ic=(DATA(m:-1:m-(Ns-1),3))’;

else %имеющееся количество выборок недостаточно для выполнения расчета ДПФ -> дополняем нулевыми значениями

Ia=[(DATA(m:-1:1,1))’,zeros(1,Ns-m)];

Ib=[(DATA(m:-1:1,2))’,zeros(1,Ns-m)];

Ic=[(DATA(m:-1:1,3))’,zeros(1,Ns-m)];

end

%%%%%%% вычисление ДПФ %%%%%%%%%%

[Ia50(m)]=fourier_calc(Ia,1);

[Ib50(m)]=fourier_calc(Ib,1);

[Ic50(m)]=fourier_calc(Ic,1);

end

%%%%%%%%% Запись в Comtrade %%%%%%%%%%

DATA_ctr=[DATA, Ia50′, Ib50′, Ic50′]; %формирование массивы выходных данных

comtrade_generator; %запуск сценария генератора comtrade

Запустим «Protection.m» на расчет и откроем осциллограмму «Osc», записанную MATLAB, в программе для просмотра и анализа осциллограмм «KIWI-Viewer». Как видно из рисунка 3 действующие значения, рассчитанные в программе MATLAB и в программе «KIWI-Viewer», совпадают.

Рисунок 3. Осциллограмма "Osc"

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

Строго говоря, приведенные формулы расчета действующего значения дают достоверный результат в случае, если основная частота сети совпадает с фундаментальным значением 50 Гц. При несовпадении данных частот, приведенный порядок расчета действующих значений не будет корректен и приведет к появлению погрешностей. В устройствах РЗА, использующих постоянную частоту дискретизации, применяют специальные мероприятия, позволяющие решить указанную проблему. Например — программную передискретизацию измеренных сигналов. А мы продолжим наше моделирования используя тестовые сигналы частотой 50 Гц.

Приложение к материалу можно скачать здесь.

Ссылка на основную публикацию
Adblock detector