sinusoida o zmiennej częstotliwości

Użytkownik Seba napisał:

A wrzuć ten skrypt/programik

To jest napisane w Matlabie

function sf=f0tof1(f0,f1,fp,tk) % sf=f0tof1(f0,f1,fp,tk) %sf - wektor sygnału o zmieniającej się częstotliwości %f0 - częstotliwość początkowa %f1 - częstotliwość końcowa %fp - częstotliwość próbkowania %tk - czas symulacji n=tk*fp t=[0:1/fp:tk-1/fp]; f=linspace(f0,f1,n); %wektor n-elementów rozłożonych liniowo od f0 do f1 sf=sin(2*pi*f.*t);

Doszedłem do wniosku że należy również modyfikować fazę sygnału ale na razie nie potrafię znaleźć w jaki sposób

Dnia 17 Apr 2005 08:46:33 +0200, Seba snipped-for-privacy@op.pl napisał:

Gwoli scislosci nie bedzie to sinusoida ;)

Zależy jak patrzeć na częstotliwość we wzorze:

y=sin(2*pi*f*t)

Jeżeli algebraicznie to jest wszystko OK, bieżesz po prostu

f=linspace(f0,f1,n)

Jak zaś patrzeć na to przez pryzmat tzw. częstotliwości chwilowej (instantaneous frequency), to świat staje się ciekawy:

- faza chwilowa w twoim wzorze: fi=(2*pi*f*t)

- częstotliwość chwilową zdefiniujmy: fi=d(fi/2/pi)/dt, co daje

fi = d(f*t)/dt = f'*t+f*t' = f'*t + f

- dla f=const, f'==0, to fc==f

- u ciebie zaś d(f)/dt == const więc dla czasu symulacji t=0->tk i f=fp->fk mamy

fc(t==0) = f'*0 + f(t==0) = fp fc(t==tk) = f'*tk + f(t==tk) = f'*tk +fk

jeżeli chcesz aby było: fck = fc(t==tk) = 40, to muszisz rozwiązać równanie (znależć fk):

fck = f'*tk + fk, gdzie f'=(fk-fp)/(tk-tp)

więc dla twoich danych fk = 30

Popatrz na wyniku skryptu

clear all

dt=1e-3; tk=1; f0=20; f1=40;

f1=(f1*tk+f0)/(tk+1)

t=0:dt:2*tk; % 2*tk, bo chcę dla t>tk f=const=f1 i0=length(t); i=length(find(t<=tk)); f=linspace(f1,f1,i0); f(1:i)=linspace(f0,f1,i);

fi=2*pi*f.*t; y=sin(fi);

% Z definicji f_d=gradient(fi)/(2*pi)/dt;

% hilbert h=hilbert(y); fi_h=unwrap(angle(h)); f_h=gradient(fi_h)/(2*pi)/dt;

% zero crossing i=find(y(1:end-1).*y(2:end)<=0); t_z=t(i)+y(i)./(y(i)-y(i+1))*dt; dt_z=diff(t_z); f_z=1./dt_z/2; t_z=t_z(1:end-1)+dt_z/2;

f_z=interp1(t_z,f_z,t); t_z=t;

% Kaiser z=y; z1=gradient(z); z2=gradient(z1); psi_y=z1.*conj(z1)-(z.*conj(z2)+z2.*conj(z))/2; z=gradient(y); z1=gradient(z); z2=gradient(z1); psi_z=z1.*conj(z1)-(z.*conj(z2)+z2.*conj(z))/2; f_k=real(asin(sqrt(psi_z./psi_y))/(2*pi)/dt);

plot( t,f, t_z,f_z, t,f_h, t,f_d, t,f_k )

Dziękuje za pomoc a szczególnie za wyjaśnienie podstaw matematycznych Dokładnie o to mi chodziło Sebastian