Быстрое преобразование Фурье с прореживанием по времени

Чувашский государственный университет имени И.Н. Ульянова

Факультет информатики и вычислительной техники

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Курсовая  работа

по дискретной математике на тему:

«Быстрое  преобразование Фурье с прореживанием по времени»

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Выполнила: студентка группы 42-09

                                                                                                             Зайкина Юлия

Проверил  преподаватель :

                                                                                                             Андреева А.А..

 

 

 

 

 

 

                             

 

 

Чебоксары 2011

Содержание:

 

Содержание: 2

Введение 3

Задание 3

1. Математическое описание  задачи 4

1.1 Определение и свойства  дискретного преобразования Фурье 4

1.2 Математическое описание  БПФ с прореживанием по времени 8

2. Разработка  алгоритмов  и программы для решения практической  задачи 12

2.1 Представление данных 12

2.2 Алгоритм быстрого  преобразования Фурье с прореживанием  по времени 12

3. Анализ результатов  моделирования алгоритма БПФ  с прореживанием по времени 14

3.1 Обращение к программе 14

3.2 Входные и выходные  параметры 14

3.3 Примеры работы программы 14

Выводы 17

Список литературы 18

Приложение 19

 

 

Введение

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

В настоящее  время во многих отраслях науки и  техники, таких как радиолокация, сейсмография, связь, радиоастрономия, медицинская электроника, цифровая обработка сигналов (ЦОС) приобретает все более широкое применение.

Важный момент применения ЦОС - это выбор разумно построенного алгоритма среди огромного разнообразия алгоритмов дискретного преобразования: быстрых алгоритмов свертки или быстрого преобразования Фурье. Общий подход к построению быстрых преобразований основан на факторизации матриц преобразования.

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

●  использование преобразований для выделения характерных признаков сигнала при спектральном анализе;

            ●  кодирование (сжатие изображений);

●  косвенное вычисление линейной свертки, позволяющее проводить ЦОС более эффективно.

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

При цифровой обработке сигнал невозможно представить непрерывной функцией, ее представляют набором значений, и используют дискретное преобразование Фурье. Алгоритмов такого преобразования достаточно много: алгоритм Кули-Тьюки, Гуда-Томаса, Герцеля и т.д. В данной работе будет рассмотрен алгоритм быстрого преобразования Фурье с прореживанием  по времени, применяемый для спектрального анализа.

Задание

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

 

1. Математическое описание задачи

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

Сигнал - это физический процесс, отображающий некоторую информацию. Сигнал может  быть математически описан некоторой функцией определенного типа.

Например:

1. одномерные сигналы описываются вещественной или комплексной функцией x_a(t), которая определена на интервале вещественной оси t'<=t<=t";

2. аналоговые сигналы - непрерывной или кусочно-непрерывной функцией x_a(t);

3. дискретные сигналы - решетчатыми функциями, или последовательностями х(nТ),

где T=const - интервал дискретизации, n =0,1,2,3,...;

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

Определение 1. Прямое преобразование Фурье - это спектр X_a(jω) аналогового сигнала x_a(t), т.е.

                                                                               (1.1)

Согласно обратному преобразованию Фурье:

                                                                       (1.2)

Определение 2. Дискретное преобразование Фурье - это пара взаимнооднозначных преобразований:

                                                           (1.3)

 

                                                        (1.4)

где Ω=2π/(NT) – основная частота преобразования.

Т.к. e^{j2π N} =W_n - поворачивающий множитель, то отсюда ДПФ и ОДПФ:

                                           (1.5)

                                         (1.6)

где x(n) - периодическая последовательность с периодом N (период – N отсчетов), т.е. x(n)=x(n+mN), m - целое число.

Определение 3. Величина называется поворачивающим множителем.

Основным свойством этих преобразований является тот факт, что из последовательности {x} получается (при прямом преобразовании) последовательность {X}, а если потом к {X} применить обратное преобразование, то снова получится исходная последовательность {x}.

Пусть имеется звуковое или какое-то иное колебание, представленное функцией x=f(t) на отрезке[0,T]. Для компьютерной обработки выполняется дискретизация: отрезок разбивается на N-1 частей и сохраняются значения функции x₀, x₁, …, x_N-1 для N точек на границах отрезков t₀=0, t₁=T/N, t₂=2T/N, …, t_n=nT/N, …, t_N-1=T.

                                      Рис.1

В результате прямого преобразования будут получены N значений X_k. Рассмотрим  теперь обратное преобразование Фурье.

                                                

X_k в общем случае величина комплексная. Разложим ее на действительную и мнимую составляющую: X_k = Re_k + i Im_k, разложим экспоненту по формуле Эйлера на синус и косинус действительного аргумента, внесем 1/N под знак суммы и перегруппируем элемент на две суммы:

Так как все x_n – действительные числа, вторая сумма равна нулю. Отсюда:

                                                   (1.7)

Поскольку при дискретизации было взято t_n=nT/N, то можно выполнить замену n=t_nN/T. В результате формула принимает следующий вид:

                                               (1.8)

Рассмотрим теперь формулу колебаний:

, где  — циклическая частота колебаний. После разложения по формуле косинуса суммы:

                                                    (1.9)

Если в этой формуле выделить компоненты, не зависящие от t и заменить  их: и , то формула запишется следующим образом:

                                                                           (1.10)

Причем про величинам Re и Im можно однозначно восстановить амплитуду и фазу исходной гармоники.

Сравним формулы (1.8) и (1.10). Видно, что сумма (1.8) представляет собой сумму из N гармонических колебаний разной частоты, фазы и амплитуды:

                                                         (1.11)

Функция называется k-гармоникой.

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

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

1. Линейность

Если X(k) и Y(k) есть ДПФ последовательностей x(nT) и y(nT) соответственно, то ДПФ последовательности ax(nT)+by(nT), где a и b – произвольные константы, равно aX(k)+bY(k).

2. Сдвиг

Пусть X(k) – ДПФ последовательности x(nT), а последовательность y(nT) получается из последовательности x(nT) путем сдвига (в случае конечной последовательности – кругового сдвига) на n₀ отсчетов: y(nT)=x(nT+n₀T). Тогда ДПФ последовательности y(nT) равно .

Аналогичный результат справедлив для сдвига коэффициентов ДПФ. Если X(k) и Y(k) есть ДПФ последовательностей x(nT) и y(nT) соответственно и Y(k)=X(k-k₀), то .

 

                                                         Рис.2

При анализе  последовательностей конечной длины  необходимо учитывать специфический  характер временного сдвига последовательности. Так, на рис.2 а изображена конечная последовательность х (n) длиной в N отсчетов. Там же крестиками изображены отсчеты эквивалентной периодической последовательности х_p(n), имеющей то же ДПФ, что и х(n). Чтобы найти ДПФ сдвинутой последовательности х(n – n₀), причем n₀ < N, следует рассмотреть сдвинутую периодическую последовательность Хр(n — n0) и в качестве эквивалентной сдвинутой конечной последовательности (имеющей ДПФ ) принять отрезок последовательности х_p(n – n₀) в интервале 0 ≤ n ≤ N — 1. Таким образом, с точки зрения ДПФ последовательность х(n – n₀) получается путем кругового сдвига элементов последовательности х(n) на n₀ отсчетов.

         3.Свойства симметрии

Если последовательность x(nT) является действительной, то ее ДПФ удовлетворяет следующим условиям симметрии:

1.2 Математическое описание БПФ с прореживанием по времени

Пусть задана последовательность x(n) , n = 0, 1, 2, … , N - 1.

                                                                  N = 2 ^t.                                                    (1.12)

Разобьем ее на две части, выделяя  четные и нечетные элементы:

                                          {x(2*n)} = {x(0),  x(2), ... , x(N–2)} и

                                 {x(2*n+1)} = {x(1),  x(3), ... , x(N–1)}.

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

Т.к. x(n) имела длину N,  X(k) имеет период N.

                                                                                                      _(1.13)

При  можно записать

                                 ;                  (1.14)

                                                                              (1.15)

Обозначим

                                ,                                     (1.16)

                              ,                                     (1.17)

тогда

                                                             _(1.18)

                             ,                (1.19)

т.е. N-точечное ДПФ вычисляется через -точечное.

Пример.     N=8.

                                                         Рис.3

Операция  “бабочки” для ДПФ с прореживанием  по времени:

         Рис.4

 

A = a + W*b

B = a – W*b.

 

4-точечное  преобразование Фурье можно вычислить  через два 2-точечных:

                                                                       Рис.5

В результате многократного прореживания входные отсчеты располагаются  в двоично-инверсном порядке. Поэтому  двоичное реверсирование является первым шагом для ДПФ с прореживанием  по времени:                                                         Таблица 1.1

Исходные номера отсчётов		Номера отсчётов после прореживания
Десятичные	Двоичные	Двоичные	Десятичные
0	000	000	0
1	001	100	4
2	010	010	2
3	011	110	6
4	100	001	1
5	101	101	5
6	110	011	3
7	111	111	7

 

Граф  алгоритма имеет следующий вид:

                          Рис.6

W_N⁰ = 1

W₄¹ = W₈².

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

Оценим  вычислительную сложность:

комплексных сложений        ;

комплексных умножений    .

Краткие выводы:

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

 

2. Разработка  алгоритмов и программы  для решения практической задачи

2.1 Представление данных

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

Входной сигнал будет одномерным массивом с длиной N. Остальные числа, такие как счётчики, будут иметь  целый тип.

В процессе вычисления происходит замещение вновь вычисленных  данных на место старых.

2.2 Алгоритм быстрого преобразования  Фурье с прореживанием по времени

На  входе дан одномерный массив чисел, т.е. вектор Х размера . На выходе необходимо получить массив чисел (спектр) размера n, т.е. преобразованный вектор Х, содержащий результат преобразования.

Эта функция БПФ реализуется на основе алгоритма, описанного ранее в п. 1.2.

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

и соответственно для ОДПФ.

 

Листинг процедуры

% перестановка элементов  согласно двоичной инверсии

for i = 1:n,            % n раз

    j = bin_invers(i - 1,n);    % номер двоично-инверсированного числа

    if j > i,            % если элемент не обработан,

        temp = x(j + 1);     % то меняем местами

        x(j + 1) = x(i);

        x(i) = temp;

    end

end;

 

% главный цикл вычислениЯ  БПФ

m = log2(n);                             % количество операций "бабочка"

for step = 1:m,

    i = 1;           % устанавливаем коэффициенты

    j = 2^(step - 1) + 1;   

    down = j;            % верхняя граница значений I

    count = 1;     % счетчик

    lim = n/2;

    k = 0;            % степень поворачивающего множителя W

    while count <= lim,          % n/2 операции "бабочка"

        xi = x(i);           % запоминаем

        xj = x(j);           % 

        Wkoef = W(2^step,k,dirFFT);   % вычисляем W    

        x(i) = xi + Wkoef*xj;    % находим новые значениЯ

        x(j) = xi - Wkoef*xj;    % выполняем "бабочку"

 

        % вычисляем  следующие i,j,k,

        if i + 1 >= down,   % если граница превышена,

            i = j + 1;       % то вычисляем следующее

            j = i + 2^(step - 1); % преобразование

            down = j;       % заново вычисляем нижнюю границу

            k = 0;           % степень W равна 0

        else            % в противном случае все значения

            i = i + 1;       % увеличиваем

            j = j + 1;       % на 1

            k = k + 1;

        end;

        count = count + 1;     % увеличиваем счетчик на 1

    end;

end;

if dirFFT ~= 0,

    X = x/n;                        % полученный результат при обратном БПФ

else

    X = x;                         % полученный результат при прямом БПФ.

End

 

3. Анализ результатов моделирования  алгоритма БПФ с прореживанием  по времени

3.1 Обращение к программе

 

Программа может быть запущена только из среды MatLab.

Для работы программы необходимо наличие  IBM PC совместимого компьютера с графическим адаптером VGA. Необходимый объем свободной оперативной памяти − не менее 35- 40 Мб. Программа запускается и работает под управлением операционной системы Windows NT, 2000, XP.

Для запуска программы в командной  строке MatLab необходимо выбрать текущий каталог программы, набрать имя программы (ptogFFT) и нажать <Enter> , либо открыть программу в отладчике и нажать клавишу <F5>  (кнопку  ) для начала ее выполнения.

 

Программа использует 5 дополнительных функций:

1. MyFFT() – вычисление БПФ с прореживанием по времени
2. Bin_invers – модуль для вычисления двоично-инверсного числа
3. W – функция для вычисления поворачивающего множителя

4)  reads_data() – считывание входных данных: N – количество отсчетов ДПФ, count – количество синусоид, входящих в сигнал, f –частоты синусоид

5)  Spectral() – проводит спектральный анализ.

3.2 Входные и выходные параметры

 

Входными данными программы  является заданная пользователем функция  (сумма синусоид), для которой строится ее график на промежутке [0,N] и график полученных результатов вычисления БПФ, которые являются выходными данными.

3.3 Примеры работы программы

 

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

Вид программы  после вызова

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

 

 

Программа после применнения написанной функции  БПФ в обратном направлении

Программа при использовании стандартной  функции БПФ

 

Выводы

 

Успешно был реализован программный  объект, реализующий быстрый алгоритм Фурье с прореживанием по времени .Приобретены навыки по Цифровой Обработке  Сигналов(ЦОС) и  реализация их в  программной среде MatLab.

 

Список  литературы

 

1. Андреева А. А. Быстрое преобразование Фурье. Конспект лекций. – Чебоксары: ЧГУ,2002. –32 с

 

2. Цифровая обработка сигналов / А. Б. Сергиенко — СПб.: Питер, 2002. — 608 с

 

 

3. Гольденберг Л. М. и др. Цифровая обработка  сигналов. –М: Радио и связь, 1985. – 312 с.

 

Приложение

Текст программы progFFT

 

function varargout = progFFT(varargin)

% PROGFFT M-file for progFFT.fig

%      PROGFFT, by itself, creates a new PROGFFT or raises the existing

%      singleton*.

%

%      H = PROGFFT returns the handle to a new PROGFFT or the handle to

%      the existing singleton*.

%

%      PROGFFT('CALLBACK',hObject,eventData,handles,...) calls the local

%      function named CALLBACK in PROGFFT.M with the given input arguments.

%

%      PROGFFT('Property','Value',...) creates a new PROGFFT or raises the

%      existing singleton*.  Starting from the left, property value pairs are

%      applied to the GUI before progFFT_OpeningFunction gets called.  An

%      unrecognized property name or invalid value makes property application

%      stop.  All inputs are passed to progFFT_OpeningFcn via varargin.

%

%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one

%      instance to run (singleton)".

%

% See also: GUIDE, GUIDATA, GUIHANDLES

 

% Edit the above text to modify the response to help progFFT

 

% Last Modified by GUIDE v2.5 08-Jun-2011 21:55:20

 

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;

gui_State = struct('gui_Name',       mfilename, ...

                   'gui_Singleton',  gui_Singleton, ...

                   'gui_OpeningFcn', @progFFT_OpeningFcn, ...

                   'gui_OutputFcn',  @progFFT_OutputFcn, ...

                   'gui_LayoutFcn',  [] , ...

                   'gui_Callback',   []);

if nargin & isstr(varargin{1})

    gui_State.gui_Callback = str2func(varargin{1});

end

 

if nargout

    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});

else

    gui_mainfcn(gui_State, varargin{:});

end

% End initialization code - DO NOT EDIT

 

 

% --- Executes just before progFFT is made visible.

function progFFT_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject    handle to figure

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

% varargin   command line arguments to progFFT (see VARARGIN)

 

% Choose default command line output for progFFT

handles.output = hObject;

 

% Update handles structure

guidata(hObject, handles);

 

% UIWAIT makes progFFT wait for user response (see UIRESUME)

% uiwait(handles.figure1);

 

 

% --- Outputs from this function are returned to the command line.

function varargout = progFFT_OutputFcn(hObject, eventdata, handles)

% varargout  cell array for returning output args (see VARARGOUT);

% hObject    handle to figure

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Get default command line output from handles structure

varargout{1} = handles.output;

 

 

% --- Executes during object creation, after setting all properties.

function edit1_CreateFcn(hObject, eventdata, handles)

% hObject    handle to edit1 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    empty - handles not created until after all CreateFcns called

 

% Hint: edit controls usually have a white background on Windows.

%       See ISPC and COMPUTER.

if ispc

    set(hObject,'BackgroundColor','white');

else

    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

 

 

function edit1_Callback(hObject, eventdata, handles)

% hObject    handle to edit1 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hints: get(hObject,'String') returns contents of edit1 as text

%        str2double(get(hObject,'String')) returns contents of edit1 as a double

 

 

% --- Executes during object creation, after setting all properties.

function edit2_CreateFcn(hObject, eventdata, handles)

% hObject    handle to edit2 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    empty - handles not created until after all CreateFcns called

 

% Hint: edit controls usually have a white background on Windows.

%       See ISPC and COMPUTER.

if ispc

    set(hObject,'BackgroundColor','white');

else

    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

 

 

 

function edit2_Callback(hObject, eventdata, handles)

% hObject    handle to edit2 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hints: get(hObject,'String') returns contents of edit2 as text

%        str2double(get(hObject,'String')) returns contents of edit2 as a double

 

 

% --- Executes during object creation, after setting all properties.

function edit3_CreateFcn(hObject, eventdata, handles)

% hObject    handle to edit3 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    empty - handles not created until after all CreateFcns called

 

% Hint: edit controls usually have a white background on Windows.

%       See ISPC and COMPUTER.

if ispc

    set(hObject,'BackgroundColor','white');

else

    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

 

 

function edit3_Callback(hObject, eventdata, handles)

% hObject    handle to edit3 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hints: get(hObject,'String') returns contents of edit3 as text

%        str2double(get(hObject,'String')) returns contents of edit3 as a double

 

 

% --- Executes on button press in pushbutton1.

function pushbutton1_Callback(hObject, eventdata, handles)

% hObject    handle to pushbutton1 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

global N count f % задание глобальных переменных

reads_data(handles);

Spectral(handles);

 

function mutual_exclude(off)

set(off,'Value',0)

 

 

% --- Executes on button press in pushbutton2.

function pushbutton2_Callback(hObject, eventdata, handles)

% hObject    handle to pushbutton2 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

close;

 

 

% --- Executes on button press in radiobutton1.

function radiobutton1_Callback(hObject, eventdata, handles)

% hObject    handle to radiobutton1 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hint: get(hObject,'Value') returns toggle state of radiobutton1

off = [handles.radiobutton2];

mutual_exclude(off)

 

 

% --- Executes on button press in radiobutton2.

function radiobutton2_Callback(hObject, eventdata, handles)

% hObject    handle to radiobutton2 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hint: get(hObject,'Value') returns toggle state of radiobutton2

off = [handles.radiobutton1];

mutual_exclude(off)

off = [handles.radiobutton4];

mutual_exclude(off)

off = [handles.radiobutton3];

mutual_exclude(off)

 

 

% --- Executes on button press in radiobutton3.

function radiobutton3_Callback(hObject, eventdata, handles)

% hObject    handle to radiobutton3 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hint: get(hObject,'Value') returns toggle state of radiobutton3

off = [handles.radiobutton4];

mutual_exclude(off)

off = [handles.radiobutton2];

mutual_exclude(off)

 

 

% --- Executes on button press in radiobutton4.

function radiobutton4_Callback(hObject, eventdata, handles)

% hObject    handle to radiobutton4 (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

 

% Hint: get(hObject,'Value') returns toggle state of radiobutton4

off = [handles.radiobutton3];

mutual_exclude(off)

off = [handles.radiobutton2];

mutual_exclude(off)

 

Текст функции myFFT

 

function X = myFFT(x,n,dirFFT);

%--------------------------------------------------------------------------

% Функция вычисления ДПФ  (Дискретное преобразование Фурье)  или ОДПФ 

% (Обратного Дискретного  Преодразования Фурье)на основе

% алгоритма БПФ (Быстрое  преобразование Фурье) с прореживанием  по времени

% (Decimation in time)

%--------------------------------------------------------------------------

% Параметры:

%   x        матрица значений входной функции

%   n        длина преобразования

%   dirFFT   направление преобразования:

%                                        =0 - прямое,

%                                       ~=0 - обратное

% Результаты:

%   а        результат преобразования 

%--------------------------------------------------------------------------

 

% перестановка элементов  согласно двоичной инверсии

for i = 1:n,            % n раз

    j = bin_invers(i - 1,n);    % номер двоично-инверсированного числа

    if j > i,            % если элемент не обработан,