Форум программистов, компьютерный форум, киберфорум
Наши страницы
Matlab
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.82/11: Рейтинг темы: голосов - 11, средняя оценка - 4.82
Kerim_Geophysic
0 / 0 / 0
Регистрация: 08.01.2016
Сообщений: 160
1

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

10.02.2016, 00:32. Просмотров 2149. Ответов 12
Метки нет (Все метки)

есть исходный сигнал 'w5'. Пусть это звуковой импульс. Известно, что поглощение спектра этого сигнала с расстоянием происходит по следующему закону:
Matlab M
1
H(f,1)=exp(-a.*f.^y.*L(1,n));
где: а-коэффициент поглощения, у-степень (обычно равна 2), f-частота спектра, L- расстояние. То есть поглощение частотно-зависимое.

я пишу код, в котором хочу разложить в спектр исходный сигнал для разных расстояний L с учетом частотно-зависимого поглощения:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
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) применялась к левой и правой стороне спектра симметрично, чтобы затем сделать обратное преобразование Фурье и проследить, как меняется сигнал с расстоянием.

есть предложения, как это реализовать?
0
Миниатюры
Прямое и обратное Фурье преобразование  
Вложения
Тип файла: rar Фурье.rar (114.6 Кб, 3 просмотров)
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
10.02.2016, 00:32
Ответы с готовыми решениями:

Двумерное прямое и обратное вейвлет-преобразование Хаара, матлаб
Подскажите пожалуйста, как сделать двумерное прямое и обратное вейвлет-преобразование Хаара в...

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

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

Обратное преобразование фурье в матлабе
Цена договорная Требуется взять обратное преобразование фурье от оределенного интеграла, нужен не...

Обратное преобразование Фурье, отсекание повторяющегося сигнала
Здравствуйте, Мой вопрос похож на тот, котрый уже тут задавался про затухание сигнала, но на тот...

12
R2D2
908 / 816 / 113
Регистрация: 23.11.2012
Сообщений: 2,418
10.02.2016, 08:48 2
Лучший ответ Сообщение было отмечено Kerim_Geophysic как решение

Решение

Kerim_Geophysic, дык просто скопируйте одну половину вектора во вторую. Только при этом лучше жестко задать размерность fft таким образом:
Matlab M
1
2
3
4
5
6
7
L = length(w);
nfft = 2^nextpow2(L);
S = fft(w, nfft);
% Потом Ваш код
% ...
% Копирование результатов
spec(:,1:nfft/2) = spec(:,nfft/2:-1:1);
1
Nick07
426 / 340 / 35
Регистрация: 17.07.2013
Сообщений: 1,808
10.02.2016, 09:02 3
В любом случае строку № 8 из цикла корректнее/надо вынести.
1
Kerim_Geophysic
0 / 0 / 0
Регистрация: 08.01.2016
Сообщений: 160
10.02.2016, 15:03  [ТС] 4
R2D2, помогла Ваша идея с копированием, я немного переделал, под свою задачу и учел тот факт, что ось симметрии находится на уровне nfft/2+1, в левой части спектра на одну гармонику оказывается больше (вроде нулевая частота, те постоянная составляющая).
вот конечный код:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
1
fft(w5)
то я получаю точно такие же числа, в том же порядке, что и в KomplSpecr(:,1), но при этом я могу высчитать обратное Фурье и получить исходный сигнал.
То есть числа одинаковые, но от своего кода я не могу сделать обратное Фурье, а от fft(w5) обратное Фурье получается хорошо. При чем если скопировать данные переменной с fft(w5) в переменную с KomplSpecr(:,1), то обратное Фурье получается.

в чем проблема?
прикрепляю файлы с кодом и рабочим листом
0
Миниатюры
Прямое и обратное Фурье преобразование  
Вложения
Тип файла: rar Фурье1.rar (85.4 Кб, 4 просмотров)
10.02.2016, 15:03
R2D2
908 / 816 / 113
Регистрация: 23.11.2012
Сообщений: 2,418
11.02.2016, 01:28 5
Лучший ответ Сообщение было отмечено Kerim_Geophysic как решение

Решение

Цитата Сообщение от Kerim_Geophysic Посмотреть сообщение
то у меня получаются комплексные числа
Лобовое решение:
y = ifft(..., 'symmetric')
А, вообще, конечно, надо делать по уму. И строго проверять симметричность. Тогда и проблем быть не должно.
И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания.
1
Kerim_Geophysic
0 / 0 / 0
Регистрация: 08.01.2016
Сообщений: 160
11.02.2016, 01:50  [ТС] 6
R2D2, спасибо огромное! Вы снова мне очень помогли))
Цитата Сообщение от R2D2 Посмотреть сообщение
И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания.
так и сделаю, ато я зациклился на этих циклах.
0
R2D2
908 / 816 / 113
Регистрация: 23.11.2012
Сообщений: 2,418
11.02.2016, 02:05 7
Kerim_Geophysic, и правильно)) Циклы в матлабе - зло. Кроме случаев, когда они неизбежны))
fft, кстати, если в нее передать матрицу - вычислит БПФ по столбцам (у Вас это каналы, как я понимаю). Это я так, на всякий случай, вдруг не знали...
1
Kerim_Geophysic
0 / 0 / 0
Регистрация: 08.01.2016
Сообщений: 160
11.02.2016, 02:12  [ТС] 8
R2D2, и снова спасибо)) потому, что я хотел писать цикл для обратного преобразования по каждому столбцу (как Вы правильно заметили - каналу), но вовремя увидел Ваш ответ)) спасибо за совет))
0
R2D2
908 / 816 / 113
Регистрация: 23.11.2012
Сообщений: 2,418
11.02.2016, 17:28 9
Kerim_Geophysic, Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее). Проверяйте, не исключено, что где то косяк. Просто хотел показать принцип. Понимаю, что код не очевиден, поэтому если с чем то не разберетесь - спрашивайте.
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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)
1
R2D2
908 / 816 / 113
Регистрация: 23.11.2012
Сообщений: 2,418
11.02.2016, 17:37 10
А на рисунке - результат обратного БПФ. Тоже пришлось использовать параметр 'symmetric'.
1
Миниатюры
Прямое и обратное Фурье преобразование  
Kerim_Geophysic
0 / 0 / 0
Регистрация: 08.01.2016
Сообщений: 160
11.02.2016, 18:00  [ТС] 11
R2D2,
Цитата Сообщение от R2D2 Посмотреть сообщение
Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее)
я называю файлы с моим именем, чтобы не перепутать то, что я написал с тем, что я где-то скачал=)) ато я создаю, потом удаляю, если мне не нравится

мне 10 строка непонятна. Я прочитал, что команда flipud(Ls) отобразит симметрично строки переменной Ls. Для чего квадратные скобки там?
0
R2D2
908 / 816 / 113
Регистрация: 23.11.2012
Сообщений: 2,418
11.02.2016, 18:31 12
Kerim_Geophysic, квадратные скобки - это так же, как мы задаем обычную матрицу. К примеру
M = [1 2 3 4; 5 6 7 8];
Но здесь не числа в явном виде, а переменные. У нас спектр зеркальный вот и функцию затухания я так же делаю зеркальной. Сначала считаю одну половинку, а потом приставляю к ней вторую но "задом на перед" (если быть точним, то "сверху вниз", т.к. функцию затухания я по столбцам записал).
1
Kerim_Geophysic
0 / 0 / 0
Регистрация: 08.01.2016
Сообщений: 160
11.02.2016, 18:52  [ТС] 13
R2D2, понял, спасибо
я вчера вот так сделал:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
ну, непрофессионально, но сойдет =)
0
11.02.2016, 18:52
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
11.02.2016, 18:52

Как построить графики с нуля по х? (дискретное и обратное преобразования Фурье)
В общем, программа строит мне графики с 1 по оси х. Как сделать, чтобы они строились от 0? Вот...

Дискретное прямое вейвлет преобразование
Здравствуйте, тема такая - Дискретное прямое вейвлет преобразование(дальше мое описание{то как я...

Обратное преобразование
Добрый день! Есть матрица (Создавалась из генератора ПСП) Произвел преобразование по Шаудеру,...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
13
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Рейтинг@Mail.ru