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

@Kerim_Geophysic · Регистрация: 08.01.2016

есть исходный сигнал 'w5'. Пусть это звуковой импульс. Известно, что поглощение спектра этого сигнала с расстоянием происходит по следующему закону:

Matlab M

1	H(f,1)=exp(-a.f.^y.L(1,n));

где: а-коэффициент поглощения, у-степень (обычно равна 2), f-частота спектра, L- расстояние. То есть поглощение частотно-зависимое.

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

Matlab M

function spec=kerim_sp(w,a,y,L)
% w - вейвлет, 1 столбец
% a - коэф поглощения
% y - степень поглощения
% L - расстояние, одна строка
S=fft(w);
for n=1:length(L);
    f=1:length(S);
    H(f,1)=exp(-a.*f.^y.*L(1,n));
    spec(:,n)=S.*H(f,1);
end

проблема в том, что спектр Фурье отображается симметрично (вроде бы это та же функция, только перевернута относительно оси У) с левой и правой стороны по оси частот.
соответственно если я код в таком виде, в каком он есть, то симметрия нарушается, и восстановить сигнал (обратное преобразование Фурье) не сделать. На рисунке это то, что у меня получается (неправильно получается

)

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

есть предложения, как это реализовать?

@R2D2 · 10.02.2016, 08:48

Kerim_Geophysic, дык просто скопируйте одну половину вектора во вторую. Только при этом лучше жестко задать размерность fft таким образом:

Matlab M

L = length(w);
nfft = 2^nextpow2(L);
S = fft(w, nfft);
% Потом Ваш код
% ...
% Копирование результатов
spec(:,1:nfft/2) = spec(:,nfft/2:-1:1);

@Nick07 · 10.02.2016, 09:02

В любом случае строку № 8 из цикла корректнее/надо вынести.

@Kerim_Geophysic · 10.02.2016, 15:03 **[ТС]**

R2D2, помогла Ваша идея с копированием, я немного переделал, под свою задачу и учел тот факт, что ось симметрии находится на уровне nfft/2+1, в левой части спектра на одну гармонику оказывается больше (вроде нулевая частота, те постоянная составляющая).
вот конечный код:

Matlab M

function spec=kerim_sp(w,a,y,L)
% w - вейвлет, 1 столбец
% a - коэф поглощения
% y - степень поглощения
% L - расстояние, одна строка
x = length(w);
nfft = 2^nextpow2(x);
S = fft(w, nfft);
for n=1:length(L);
    f=1:length(S)
    H(f,1)=exp(-a.*f.^y.*L(1,n));
    spec(:,n)=S.*H(f,1);
end
spec(nfft/2+2:nfft,:) = spec(nfft/2:-1:2,:); % Копирование результатов
end

Nick07, учел Ваш совет, спасибо

возникает другая проблема, не могу сделать обратное преобразование.
тут все очень странно. Если я беру спектр по своему коду для первого столбца ( то есть спектр исходного сигнала, где L=0 и нет затухания) и пытаюсь сделать обратное Фурье преобразование, то у меня получаются комплексные числа.

Matlab M

1	ifft(KomplSpecr(:,1))

Если я буре спектр от исходного сигнала

Matlab M

fft(w5)

то я получаю точно такие же числа, в том же порядке, что и в KomplSpecr(:,1), но при этом я могу высчитать обратное Фурье и получить исходный сигнал.
То есть числа одинаковые, но от своего кода я не могу сделать обратное Фурье, а от fft(w5) обратное Фурье получается хорошо. При чем если скопировать данные переменной с fft(w5) в переменную с KomplSpecr(:,1), то обратное Фурье получается.

в чем проблема?
прикрепляю файлы с кодом и рабочим листом

@R2D2 · 11.02.2016, 01:28

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

то у меня получаются комплексные числа

Лобовое решение:
y = ifft(..., 'symmetric')
А, вообще, конечно, надо делать по уму. И строго проверять симметричность. Тогда и проблем быть не должно.
И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания.

@Kerim_Geophysic · 11.02.2016, 01:50 **[ТС]**

R2D2, спасибо огромное! Вы снова мне очень помогли))

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

И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания.

так и сделаю, ато я зациклился на этих циклах.

@R2D2 · 11.02.2016, 02:05

Kerim_Geophysic, и правильно)) Циклы в матлабе - зло. Кроме случаев, когда они неизбежны))
fft, кстати, если в нее передать матрицу - вычислит БПФ по столбцам (у Вас это каналы, как я понимаю). Это я так, на всякий случай, вдруг не знали...

@Kerim_Geophysic · 11.02.2016, 02:12 **[ТС]**

R2D2, и снова спасибо)) потому, что я хотел писать цикл для обратного преобразования по каждому столбцу (как Вы правильно заметили - каналу), но вовремя увидел Ваш ответ)) спасибо за совет))

@R2D2 · 11.02.2016, 17:28

Kerim_Geophysic, Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее). Проверяйте, не исключено, что где то косяк. Просто хотел показать принцип. Понимаю, что код не очевиден, поэтому если с чем то не разберетесь - спрашивайте.

Matlab M

clc; close all; clear;
load spectrums
spec = kerim_sp(w5,a,2,L);
a = 1e-5;
n = 2;
nfft = 2^nextpow2(length(w5));
% nfft = length(w5);
Dist = 0:100;
Ls = exp(-a*(1:nfft/2)'.^n*Dist);
Ls = [flipud(Ls); Ls];
S = fft(w5, nfft);
spec1 = (S*ones(1, length(Dist))).*Ls;
%%
figure('Units', 'Normalized', 'Position', [.05 .05 .5 .84])
surf(abs(spec1));
shading interp
colormap summer
camlight(45, 45)
view(-60, 30)

@R2D2 · 11.02.2016, 17:37

А на рисунке - результат обратного БПФ. Тоже пришлось использовать параметр 'symmetric'.

@Kerim_Geophysic · 11.02.2016, 18:00 **[ТС]**

R2D2,

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

Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее)

я называю файлы с моим именем, чтобы не перепутать то, что я написал с тем, что я где-то скачал=)) ато я создаю, потом удаляю, если мне не нравится

мне 10 строка непонятна. Я прочитал, что команда flipud(Ls) отобразит симметрично строки переменной Ls. Для чего квадратные скобки там?

@R2D2 · 11.02.2016, 18:31

Kerim_Geophysic, квадратные скобки - это так же, как мы задаем обычную матрицу. К примеру
M = [1 2 3 4; 5 6 7 8];
Но здесь не числа в явном виде, а переменные. У нас спектр зеркальный вот и функцию затухания я так же делаю зеркальной. Сначала считаю одну половинку, а потом приставляю к ней вторую но "задом на перед" (если быть точним, то "сверху вниз", т.к. функцию затухания я по столбцам записал).

@Kerim_Geophysic · 11.02.2016, 18:52 **[ТС]**

R2D2, понял, спасибо
я вчера вот так сделал:

Matlab M

function spec=kerim_sp_atten(w,a,y,L)
% w - вейвлет, 1 столбец
% a - коэф поглощения
% y - степень поглощения
% L - расстояние, одна строка
x = length(w);
nfft = 2^nextpow2(x);
S = fft(w, nfft);
k=repmat(L,length(S),1);
f=1:length(S)
g=repmat(f',1,length(L));
H=exp(-a.*g.^y.*k);
spec=repmat(S,1,length(L)).*H;
spec(nfft/2+2:nfft,:) = spec(nfft/2:-1:2,:); % Копирование результатов
end

ну, непрофессионально, но сойдет =)

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

Решение

Решение

@Nick07 529 / 434 / 46 Регистрация: 17.07.2013 Сообщений: 2,223
	10.02.2016, 09:02	3
	В любом случае строку № 8 из цикла корректнее/надо вынести. 1

@R2D2 938 / 838 / 121 Регистрация: 23.11.2012 Сообщений: 2,489
	11.02.2016, 01:28	5
	Сообщение было отмечено Kerim_Geophysic как решение Решение Сообщение от Kerim_Geophysic то у меня получаются комплексные числа Лобовое решение: y = ifft(..., 'symmetric') А, вообще, конечно, надо делать по уму. И строго проверять симметричность. Тогда и проблем быть не должно. И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания. 1

@Kerim_Geophysic 2 / 2 / 0 Регистрация: 08.01.2016 Сообщений: 483
	11.02.2016, 01:50 [ТС]	6
	R2D2, спасибо огромное! Вы снова мне очень помогли)) Сообщение от R2D2 И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания. так и сделаю, ато я зациклился на этих циклах. 0

@R2D2 938 / 838 / 121 Регистрация: 23.11.2012 Сообщений: 2,489
	11.02.2016, 02:05	7
	Kerim_Geophysic, и правильно)) Циклы в матлабе - зло. Кроме случаев, когда они неизбежны)) fft, кстати, если в нее передать матрицу - вычислит БПФ по столбцам (у Вас это каналы, как я понимаю). Это я так, на всякий случай, вдруг не знали... 1

@Kerim_Geophysic 2 / 2 / 0 Регистрация: 08.01.2016 Сообщений: 483
	11.02.2016, 02:12 [ТС]	8
	R2D2, и снова спасибо)) потому, что я хотел писать цикл для обратного преобразования по каждому столбцу (как Вы правильно заметили - каналу), но вовремя увидел Ваш ответ)) спасибо за совет)) 0

@R2D2 938 / 838 / 121 Регистрация: 23.11.2012 Сообщений: 2,489
	11.02.2016, 17:37	10
	А на рисунке - результат обратного БПФ. Тоже пришлось использовать параметр 'symmetric'. Миниатюры 1

@Kerim_Geophysic 2 / 2 / 0 Регистрация: 08.01.2016 Сообщений: 483
	11.02.2016, 18:00 [ТС]	11
	R2D2, Сообщение от R2D2 Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее) я называю файлы с моим именем, чтобы не перепутать то, что я написал с тем, что я где-то скачал=)) ато я создаю, потом удаляю, если мне не нравится мне 10 строка непонятна. Я прочитал, что команда flipud(Ls) отобразит симметрично строки переменной Ls. Для чего квадратные скобки там? 0

@R2D2 938 / 838 / 121 Регистрация: 23.11.2012 Сообщений: 2,489
	11.02.2016, 18:31	12
	Kerim_Geophysic, квадратные скобки - это так же, как мы задаем обычную матрицу. К примеру M = [1 2 3 4; 5 6 7 8]; Но здесь не числа в явном виде, а переменные. У нас спектр зеркальный вот и функцию затухания я так же делаю зеркальной. Сначала считаю одну половинку, а потом приставляю к ней вторую но "задом на перед" (если быть точним, то "сверху вниз", т.к. функцию затухания я по столбцам записал). 1