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

[_20_]

разделы "%Вот здесь я пытаюсь изменить вектор времени." и "% Вот здесь я дополнительно избегаю зеркального эффекта." отличаются использованием в трансформации F1 вместо F.
F - это непосредственно коэфиценты после трансформации, F1 - это исправленные коэфиценты:

F = fft(S);                                                        % Это фурье образ, здесь всё в порядке. F - это коэфиценты разложения.
F2 = abs(F/size(S,2));                                      %Нормированный модуль коэфицентов (делёный на количество значений в сигнале (длина вектора S))
F1 = F2(1:size(S,2)/2+1);                               % Берём половину от нормированного результата трансформации (полученого строчкой выше)
F1(2:end-1) = 2*F1(2:end-1);                         % Из этой половины все, за исключением первого значения, удваиваем.

-- 26.09.2018, 00:42 --

"Signal after transformation" - это результат обратной трансформации в которой я использовал те же самые значения времени, что и при прямой трансформации. В этом случае даже было не обязательно коректировать значения прямой трансформации и бороться с зеркальным эффектом (я сипользовал значения F, а не F1) и всё нормально сложилось в первоначальный сигнал.

[Pavia]

_20_ в сообщении #1341232 писал(а):

temp3(n) = temp3(n) + abs(F(k)) * exp(1i * angle(F(k))) * exp((2 * pi * 1i * (n - 1) * (k - 1))/(size(t3,2)));

Здесь надо делить не на размер массива, а на частоту дескритизации. Она в 1.5 раза меньше.

Плюс у вас растекание спектра.
Поэтому перед fft(S) надо убрать постоянную составляющую.
S = S - mean(S)

[_20_]

Pavia

Я Ваших рекомендаций не понял, но всё - таки им последовал:

S = 0.7*sin(2*pi*3*t) + sin(2*pi*6*t);
S = S - mean(S);

temp3(n) = temp3(n) + abs(F1(k)) * exp(1i * angle(F1(k))) * exp((2 * pi * 1i * (n - 1) * (k - 1))/(size(t3,2)/1.5));

Не помогло:

[arseniiv]

Кстати, зачем делать вот это: abs(F1(k)) * exp(1i * angle(F1(k)))? Это же F1(k) и есть.

-- Ср сен 26, 2018 17:14:51 --

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

О Диэдр, теорию никто не изучает и справку никто не читает.

[wrest]

_20_
Я попробовал пофурьерить сам
Подозреваю, что проблема вот тут
abs(F1(k)) * exp(1i * angle(F1(k))) * exp((2 * pi * 1i * (n - 1) * (k - 1))/(size(t3,2)/1.5));
Я делал "вещественное ДПФ" с вашим сигналом S = 0.7*sin(2*pi*3*t) + sin(2*pi*6*t);
Прямое ДПФ
$ReF(k)=\sum \limits_{i=0}^{N-1}S_i\cos(2\pi ki/N)$
$ImF(k)=-\sum \limits_{i=0}^{N-1}S_i\sin(2\pi ki/N)$
Тут соответственно $N=150;S_i=S(i\cdot \Delta t); \Delta t = 0.01$
Нормализация:
$ReF(1)=2\cdot ReF(1);ReF(N/2)=2\cdot ReF(N/2)$
$ImF(1)=2\cdot ImF(1);ImF(N/2)=2\cdot ImF(N/2)$
Обратное ДПФ (т.е. -- разложение в ряд):
$f(t)=\dfrac{1}{N}\left(\sum \limits_{i=0}^{N/2}ReF(i)\cdot \cos(2\pi i \dfrac{t}{\Delta t N})-\sum \limits_{i=0}^{N/2}ReF(i)\cdot \sin (2\pi i \dfrac{t}{\Delta t N})\right)$
Это то что вы хотели -- $f$ это функция от непрерывного времени $t$

$S(t), t \in [0;1]$:

$f(t), t \in [0;1]$:

$S(t)-f(t), t \in [0;1]$

Там мелковато, размах ошибки: $(-0.06;-0.02)$

[arseniiv]

У ТС там вроде проблема в том, что если отрезать вторую половину списка амплитуд, будет не то, и если её не отрезать, но не трактовать по-правильному как амплитуды отрицательных частот — тоже будет не то (хотя на отсчётах значения, конечно, совпадут — и то, и то понятно почему). Советы Pavia я вообще не понял, особенно множитель 1,5. У всех операций есть смысл, их нельзя делать как попало.

[wrest]

arseniiv в сообщении #1341650 писал(а):

Советы Pavia я вообще не понял, особенно множитель 1,5.

Ну интервал $1,5$ секунды, отсчеты берем через $0,01$ секунд, частота дискретизации $100$ герц, всего $150$ отсчетов. Для перехода от отсчетов к непрерывному времени надо время поделить на полтора ($\Delta t \cdot N=0,01 \cdot 150 = 1,5$)

[arseniiv]

Ладно, снимаю возражение как вызванное трудностями в прочтении постов.

[wrest]

wrest в сообщении #1341646 писал(а):

Обратное ДПФ (т.е. -- разложение в ряд):
$f(t)=\dfrac{1}{N}\left(\sum \limits_{i=0}^{N/2}ReF(i)\cdot \cos(2\pi i \dfrac{t}{\Delta t N})-\sum \limits_{i=0}^{N/2}ReF(i)\cdot \sin (2\pi i \dfrac{t}{\Delta t N})\right)$
Это то что вы хотели -- $f$ это функция от непрерывного времени $t$

Читать так:
Обратное ДПФ (т.е. -- разложение в ряд):
$f(t)=\dfrac{1}{N}\left(\sum \limits_{i=0}^{N/2}ReF(i)\cdot \cos(2\pi i \dfrac{t}{\Delta t N})-\sum \limits_{i=0}^{N/2}\color{blue}Im\color{black} F(i)\cdot \sin (2\pi i \dfrac{t}{\Delta t N})\right)$
Это то что вы хотели -- $f$ это функция от непрерывного времени $t$
Да, кстати. Здесь $i$ это индекс, а не мнимая единица, есличо :)

Красный цвет зарезервирован под модераторскую правку. Исправлено на синий. // Lia

[_20_]

Здравствуйте снова.
Я взял не большую паузу, но никак не могу понять, в чём ошибка. Вроде делаю всё правильно, проверьте?
Не получается обратное фурье преобразование даже на том же самом отрезке времени:

step = 0.01; %in seconds
SignalLength = 1.5; %in seconds
t = (0:step:SignalLength-step);
S = 0.7*sin(2*pi*3*t) + sin(2*pi*6*t);
a = zeros(1,size(t,2));
b = a;
F = complex(a,0);
for w = 1:size(t,2)
    for j = 1:size(t,2)
        a(w) = a(w) + S(j) * cos(2*pi*(w-1)*(j-1)/size(t,2));
        b(w) = b(w) - S(j) * sin(2*pi*(w-1)*(j-1)/size(t,2));
        F(w) = a(w) + 1i * b(w);
    end
end
F(1) =  2*F(1);
F(size(t,2)/2) = 2 * F(size(t,2)/2);
f = zeros(size(S,2));
for k = 1:size(t,2)
    for w = 1:size(t,2)/2   
        f(k) = f(k) + 1/size(t,2) * (real(F(w)) * cos(2*pi*(w-1)*(k-1)/size(t,2)/step)...
            - imag(F(w)) * sin(2*pi*(w-1)*(k-1)/size(t,2)/step));
    end
end

[arseniiv]

Вы бы имена в коде подробнее документировали, что ли, или хотя бы сказали, что с чем не сходится, ну и проверили соответствие формулам. И например для не разбирающихся в матлабе что делает S(j)? Если понимаю код правильно, проблема могла бы быть во множителях внутри синуса и косинуса.

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

Ещё строку F(w) = a(w) + 1i * b(w); следует вынести во внешний цикл for w, потому что лишние присвоения во внутреннем не дают никакого эффекта.

[wrest]

_20_ в сообщении #1342535 писал(а):

Вроде делаю всё правильно, проверьте?

Я не понимаю матлаб
Тоже не понял как работает S(j)
Что такое size(t,2)? Это целое, N? Ну и назвади бы тогда N=size(t,2)...
Как работает конструкция с двумя делениями подряд 2*pi*(w-1)*(k-1)/size(t,2)/step ? И зачем оно там, если вы вычисляете $f_k$ а не $f(t)$ -- похоже тут и ошибка.
Вы делаете циклы от 1, но внутри от индекса отнимаете единицу. Матлаб не поддержтвает индексы от нуля?
Замените real(F(w)) и imag(F(w)) на a(w) и b(w) соотвественно, у вас же это вычисляется ранее. И следите за индексами (или делайте от нуля или не забывайте отнимать единицу).
То есть: сделайте точно по формулам вещественного Фурье что я вам написал. Потом будете дорабатывать и анализировать почему комплексное Фурье не получается.

[_20_]

Спасибо, что не бросаете тему, я сейчас исправлюсь.
wrest Вроде делаю всё именно по тем формулам, что Вы написали, но, видимо, делаю что - то неправильно и не могу понять, что именно.
Да, в матлабе принята индексация масивов от еденицы, как в фортране.
Конструкция с двумя делениями работает так же, как конструкция с двумя умножениями, сначало выполняется деление на первое число, а потом результат делится на второе число.

Перевод с матлаба на человеческий:

step = 0.01; %in seconds
Это я принимаю шаг приращения времени для вычисления входящего сигнала, наверное правильно назвать дискретизация. Оно в секундах, мне так будет удобнее.

SignalLength = 1.5; %in seconds
Это длинна входящего сигнала, тоже в секундах.

t = (0:step:SignalLength-step);
Это создание значений времени. Является одномерным масивом, в матлабе такое принято называть вектором, но я постараюсь здесь называть именно одномерный массив.

S = 0.7*sin(2*pi*3*t) + sin(2*pi*6*t);
Это создание массива интенсивности входящего сигнала. Тоже одномерный массив.

a = zeros(1,size(t,2));
Создание пустого одномерного массива длинной в количество точек входящего сигнала. Будет использоваться для содержания действительной части фурье образа.

b = a;
Создание пустого одномерного массива длинной в количество точек входящего сигнала. Будет использоваться для содержания мнимой части фурье образа.

F = complex(a,0);
Создание пустого одномерного массива комплексных чисел длинной в количество точек входящего сигнала. Будет использоваться для содержания всего фурье образа.

for w = 1:size(t,2)
Цикл. w будет пробегать все значения от одного до количества пунктов одномерного массива времени.

for j = 1:size(t,2)
Цикл. j будет пробегать все значения от одного до количества пунктов одномерного массива времени.

a(w) = a(w) + S(j) * cos(2*pi*(w-1)*(j-1)/size(t,2));
Вычисление действительной части фурье образа. S(j) - это j-ый элемент одномерного массива S, входящего сигнала.

b(w) = b(w) - S(j) * sin(2*pi*(w-1)*(j-1)/size(t,2));
Вычисление мнимой части фурье образа. S(j) - это j-ый элемент одномерного массива S, входящего сигнала.

F(w) = a(w) + 1i * b(w);
Вычисление фурье образа.

end
Конец внутреннего цикла.

end
Конец внешнего цикла.

F(1) = 2*F(1);
Удвоение первого элемента фурье образа.

F(size(t,2)/2) = 2 * F(size(t,2)/2);
Удвоение среднего элемента фурье образа (то есть того элемента, который находится по середине одномерного массива фурье образа).

f = zeros(size(S,2));
Создание пустого одномерного массива длинной в количество точек входящего сигнала. Будет использоваться для содержания сигнала, полученного обратной трансформацией.

for k = 1:size(t,2)
Цикл. k будет пробегать все значения от одного до количества пунктов одномерного массива времени.

for w = 1:size(t,2)/2
Цикл. w будет пробегать все значения от одного до половины от количества пунктов одномерного массива времени.

f(k) = f(k) + 1/size(t,2) * (real(F(w)) * cos(2*pi*(w-1)*(k-1)/size(t,2)/step)...
- imag(F(w)) * sin(2*pi*(w-1)*(k-1)/size(t,2)/step));
Вычисление обратной трансформации фурье. k - это типа время, w - частота. То есть для каждого момента времени я складываю по всем частотам с их коэфицентами.

end
Конец внутреннего цикла.

end
Конец внешнего цикла.

[wrest]

_20_ в сообщении #1343192 писал(а):

a = zeros(1,size(t,2));
Создание пустого одномерного массива длинной в количество точек входящего сигнала.

Что значит "2" в size(t,2)?

_20_ в сообщении #1343192 писал(а):

f(k) = f(k) + 1/size(t,2) * (real(F(w)) * cos(2*pi*(w-1)*(k-1)/size(t,2)/step)...
- imag(F(w)) * sin(2*pi*(w-1)*(k-1)/size(t,2)/step));
Вычисление обратной трансформации фурье. k - это типа время, w - частота.

Что значит "типа время"?
Вот же формулы для непрерывного времени:

wrest в сообщении #1341684 писал(а):

Обратное ДПФ (т.е. -- разложение в ряд):
$f(t)=\dfrac{1}{N}\left(\sum \limits_{i=0}^{N/2}ReF(i)\cdot \cos(2\pi i \dfrac{t}{\Delta t N})-\sum \limits_{i=0}^{N/2}\color{blue}Im\color{black} F(i)\cdot \sin (2\pi i \dfrac{t}{\Delta t N})\right)$
Это то что вы хотели -- $f$ это функция от непрерывного времени $t$

Где вы там видите внешний цикл? Аргумент у синусов\косинусов или
$2\pi i \cdot \dfrac{t}{\Delta t} \cdot \dfrac{1}{N}$ - если вы хотите считать для любого $t$
Или $2\pi i \cdot k \cdot \dfrac{1}{N}$ - если вы восстанавливаете ваши 0-й, 1-й, ... 149-й отсчет и уходите от времени... То есть $k$-й отсчет соответствует моменту времени $t=k \Delta t$
Вы хотели восстанавливать значения функции которые попали между отсчетами при дискретизации, например между 1 и 2 отсчетом, тогда у вас получается отсчет за номером 1,5 к примеру, вот для этого вы в аргументы синусов\косинусов и подставляете вместо номера отсчета $k$ значение времени $t$, в виде $k=\dfrac{t}{\Delta t}$

[_20_]

Спасибо за ответы. У меня возникли некоторые сложности, но, вроде я сам разобрался. Осталась пара вопросов.

wrest
По Вашим формулам обратное ДПФ даёт амплитуду сигнала в два раза меньше. Надо или на два умножать, ну или делить не на N, а на N/2.

Метод хорошо работает не для всех длинн входящего сигнала. Вот пример:
https://ibb.co/nwtvEU
Верхний ряд - время t проходин от 0 до 1,

SignalLength = 1; 
t = (0:step:SignalLength-step);
S = 0.7*sin(2*pi*3*t) + sin(2*pi*6*t);

нижний ряд - время проходит от 0 до 1,1

SignalLength = 1.1;
t = (0:step:SignalLength-step);
S = 0.7*sin(2*pi*3*t) + sin(2*pi*6*t);

Как видно, появляются постоянные составляющие сигнала. Отчего это может быть и как от них избавиться?

Второй вопрос, что делать, если длина сигнала нечётная?

Вы, при нормализации, удваиваете:

$ ReF(N/2) = 2 ReF(N/2)$

Что делать, если N/2 не является целым числом?