Функция вейвлет-преобразования Морле возвращает бессмысленный график

Я написал функцию matlab (версия 7.10.0.499 (R2010a)) для оценки входящего сигнала FT и вычисления вейвлета Морле для сигнала. У меня есть похожая программа, но мне нужно было сделать ее более читаемой и ближе к математическому жаргону. Выходной график должен быть 2D-график с цветом, показывающим интенсивность частот. Мой сюжет, кажется, все частоты одинаковы в одно и то же время. Программа действительно делает fft в строке времени для каждой частоты, поэтому я предполагаю, что другой способ взглянуть на это заключается в том, что одна и та же линия повторяется на каждом шаге в моем цикле for. Проблема в том, что я проверил исходную программу, которая возвращает правильный сюжет, и я не могу найти никакой разницы, кроме того, что я назвал значения и как я организовал код.

function[msg] = mile01_wlt(FT_y, f_mn, f_mx, K, N, F_s)
%{
Fucntion to perform a full wlt of a morlet wavelett.
optimization of the number of frequencies to be included. 
FT_y satisfies the FT(x) of 1 envelope and is our ft signal.
f min and max enter into the analysis and are decided from 
the f-image for optimal values.
While performing the transformation there are different scalings
on the resulting "intensity".
Plot is made with a 2D array and a colour code for intensity. 
version 05.05.2016
%}

%--------------------------------------------------------------%
%{
tableofcontents:
    1: determining nr. of analysis f, prints and readies f's to be used.
    2: ensuring correct orientation of FT_y
    3:defining arrays
    4: declaring waveletdiagram and storage of frequencies
    5: for-loop over all frequencies:
    6: reducing file to manageable size by truncating time.
    7: marking plot to highlight ("randproblemer")
    8: plotting waveletdiagram
%}

%--------------------------------------------------------------%
%1: determining nr. of analysis f, prints and readies f's to be used.
    DF = floor( log(f_mx/f_mn) / log(1+( 1/(8*K) ) ) ) + 1;% f-spectre analysed
    nr_f_analysed = DF              %output to commandline
    f_step = (f_mx/f_mn)^(1/(DF-1)); % multiplicative step for new f_a
    f_a = f_mn; %[Hz] frequency of analysis
    T = N/F_s; %[s] total time sampled
    C = 2.0; % factor to scale Psi

%--------------------------------------------------------------%
%2: ensuring correct orientation of FT_y
    siz = size(FT_y);
    if (siz(2)>siz(1)) 
        FT_y = transpose(FT_y);
    end;

%--------------------------------------------------------------%    
%3:defining arrays
    t = linspace(0, T*(N-1)/N, N); %[s] timespan
    f = linspace(0, F_s*(N-1)/N, N); %[Hz] f-specter

%--------------------------------------------------------------%
%4: declaring waveletdiagram and storage of frequencies
    WLd = zeros(DF,N); % matrix of DF rows and N collumns for storing our wlt
    f_store = zeros(1,DF); % horizontal array for storing DF frequencies

%--------------------------------------------------------------%
%5: for-loop over all frequencies:
    for jj = 1:DF
        o = (K/f_a)*(K/f_a); %factor sigma
        Psi = exp(- 0*(f-f_a).*(f-f_a)); % FT(psi) for 1 envelope
        Psi = Psi - exp(-K*K)*exp(- o*(f.*f)); % correctional element
        Psi = C*Psi; %factor. not set in stone

        %next step fits 1 row in the WLd (3 alternatives)
        %WLd(jj,:) = abs(ifft(Psi.*transpose(FT_y))); 
        WLd(jj,:) = sqrt(abs(ifft(Psi.*transpose(FT_y))));
        %WLd(jj,:) = sqrt(abs(ifft(Psi.*FT_y))); %for different array sizes
                                                %and emphasizes weaker parts.
        %prep for next round
        f_store (jj) = f_a; % storing used frequencies
        f_a = f_a*f_step; % determines the next step
    end;

%--------------------------------------------------------------%
%6: reducing file to manageable size by truncating time.
    P = floor( (K*F_s) / (24*f_mx) );%24 not set in stone
    using_every_P_point = P %printout to cmdline for monitoring
    N_P = floor(N/P);
    points_in_time = N_P %printout to cmdline for monitoring
    % truncating WLd and time
    WLd2 = zeros(DF,N_P);
    for jj = 1:DF
        for ii = 1:N_P
            WLd2(jj,ii) = WLd(jj,ii*P);
        end
    end
    t_P = zeros(1,N_P);
    for ii = 1:N_P % set outside the initial loop to reduce redundancy
        t_P(ii) = t(ii*P);
    end

%--------------------------------------------------------------%
%7: marking plot to highlight boundary value problems
    maxval = max(WLd2);%setting an intensity
    mxv = max(maxval);
    % marks in wl matrix
    for jj= 1:DF
        m = floor( K*F_s / (P*pi*f_store(jj)) ); %finding edges of envelope
        WLd2(jj,m) = mxv/2; % lower limit
        WLd2(jj,N_P-m) = mxv/2;% upper limit
    end

%--------------------------------------------------------------%
%8: plotting waveletdiagram
    figure;
    imagesc(t_P, log10(f_store), WLd2, 'Ydata', [1 size(WLd2,1)]);
    set(gca, 'Ydir', 'normal');
    xlabel('Time [s]');
    ylabel('log10(frequency [Hz])');
    %title('wavelet power spectrum'); % for non-sqrt inensities
    title('sqrt(wavelet power spectrum)'); %when calculating using sqrt
    colorbar('location', 'southoutside');
    msg = 'done.';

Нет никакого сообщения об ошибке, поэтому я не уверен, что именно я делаю неправильно.
Надеюсь, я следовал всем рекомендациям. В противном случае, я прошу прощения.

редактировать:
моя вызывающая программа:
% установление параметров
N = 2^(16); % / количество точек для выборки
F_s = 3.2e6; % частоты Hz | samplings
T_t = N / F_s; % S / длина в секундах времени выборки
f_c = 2.0e5; % Hz / частота несущей волны
f_m = 8./T_t; % Гц / частота модулирующей волны
w_c = 2 * pi * f_c; % Гц / угловая частота («омега») несущей волны
w_m = 2 * pi * f_m; % Гц / угловая частота («омега») модулирующей волны

% establishing parameter arrays
t = linspace(0, T_t, N);

% function variables
T_h = 2*f_m.*t; % dimless | 1/2 of the period for square signal

% combined carry and modulated wave
% y(t) eq. 1):
y_t = 0.5.*cos(w_c.*t).*(1+cos(w_m.*t)); 
% y(t) eq. 2):
%     y_t = 0.5.*cos(w_c.*t)+0.25*cos((w_c+w_m).*t)+0.25*cos((w_c-w_m).*t);
%square wave
sq_t = cos(w_c.*t).*(1 - mod(floor(t./T_h), 2)); % sq(t)

% the following can be exchanged between sq(t) and y(t) 

plot(t, y_t)
% plot(t, sq_t)
xlabel('time [s]');
ylabel('signal amplitude');
title('plot of harmonically modulated signal with carrying wave');
% title('plot of square modulated signal with carrying wave');
figure()
hold on

% Fourier transform and plot of freq-image
FT_y = mile01_fftplot(y_t, N, F_s);
% FT_sq = mile01_fftplot(sq_t, N, F_s);

% Morlet wavelet transform and plot of WLdiagram
%determining K, check t-image
K_h = 57*4; % approximation based on 1/4 of an envelope, harmonious
%determining f min and max, from f-image
f_m = 1.995e5; % minimum frequency. chosen to showcase all relevant f
f_M = 2.005e5; % maximum frequency. chosen to showcase all relevant f
%calling wlt function.
name = 'mile'
msg = mile01_wlt(FT_y, f_m, f_M, K_h, N, F_s)
siz = size(FT_y);
if (siz(2)>siz(1)) 
    FT_y = transpose(FT_y);
end;
name = 'arnt'
msg = arnt_wltransf(FT_y, f_m, f_M, K_h, N, F_s)

Временное изображение имеет постоянную частоту, но амплитуда колеблется, напоминая гауссову кривую. Мой код возвращает резко сегментированное изображение с течением времени, где каждая точка во времени содержит только 1 частоту. Он должен отражать изменение интенсивности во всех спектрах с течением времени.
надеюсь, что поможет и спасибо!

1 ответ

  1. Я нашел ошибку. Существует 0, а не o в первом экземпляре Psi. Думаю, что я, возможно, переименую значение как sig или что-то в этом роде. кроме того, код работает. извините за неприятности там