La fonction de transformation en ondelettes de Morlet renvoie un tracé non sensuel

Question

J'ai écrit une fonction matlab (Version 7.10.0.499 (R2010a)) pour évaluer le signal FT entrant et calculer l'ondelette Morlet pour le signal. J'ai un programme similaire, mais je devais le rendre plus lisible et plus proche du jargon mathématique. Le tracé de sortie est supposé être un tracé 2D avec une couleur indiquant l'intensité des fréquences. Mon intrigue semble avoir toutes les fréquences les mêmes par moment. Le programme fait un fft par ligne de temps pour chaque fréquence, donc je suppose qu'une autre façon de le regarder est que la même ligne se répète par pas dans ma boucle for. Le problème est que j'ai vérifié avec le programme original, qui retourne l'intrigue correcte, et je ne peux pas trouver de différence au-delà de ce que j'ai nommé les valeurs et comment j'ai organisé le code.

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.'; 

Il n'y a pas de message d'erreur, donc je ne sais pas exactement ce que je fais mal. J'espère avoir suivi toutes les directives. Sinon, je m'excuse.

modifier: mon programme appelant: % paramètres d'établissement N = 2^(16); % | nombre de points à échantillonner F_s = 3,2e6; % Hz | échantillonnages fréquence T_t = N/F_s; % s | longueur en secondes du temps d'échantillonnage f_c = 2,0e5; % Hz | fréquence d'onde porteuse f_m = 8./T_t; % Hz | fréquence d'onde modulante w_c = 2 * pi * f_c; % Hz | fréquence angulaire ("oméga") de l'onde portante w_m = 2 * pi * f_m; % Hz | la fréquence angulaire (« oméga ») de modulation onde

% 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) 

L'image de temps a une fréquence constante, mais l'amplitude oscille resempling une courbe gaussienne. Mon code renvoie une image fortement segmentée dans le temps, où chaque point de temps ne contient qu'une seule fréquence. Il devrait refléter un changement d'intensité à travers les spectres au fil du temps. espérons que cela aide et merci!

Answer 1

J'ai trouvé l'erreur. Il y a un 0 plutôt qu'un o dans la première instance de Psi. Penser je vais peut-être renommer la valeur en tant que sig ou quelque chose. en outre, le code fonctionne. désolé pour le problème là