how to generate a Sine wave with changable frequency in mfile?

8 views (last 30 days)
Dear all,
I would like to generate something similler to a PLL sine signal locked on a measured sine signal. for that i have detected the Zeros of my measured, lets say, voltage signal and calculated the freuqncy from counting the number of points between two zero crossing. (the frequency of the measured signal is not constant and changing with time). by using the following equation that was achieved:
f = 1./(2*diff(Zeros));% where Zeros are my detected zerocorssing times in s.
so by knowing the following:
1- I know the frequency of my measured signal (f vector length is not equal to the length of measured signal)
2- I know my sample rate of my measured signal and my time.
3- sin(2*pi*f*t) will generate my locked sine signal to the measured sine wave signal (Voltage)
my question: why am i not getting a correct sine wave signal with the same frequency of the measured signal?
note: Zeros is a vector containing the interpolated time values when my sine wave signal corss the zero value. (so i only have frequency values at the diff(Zeros) values and not at the same sample rate of the measured signal).
Can anyone help me to condition my calculated frequency signal to have the same lenght as the measured signal (and the correct time steps) in order to make the sin(2*pi*f*t) works correctly?
Best Regards
  2 Comments
Aubai
Aubai on 11 Aug 2014
thx for trying to help my code looks like the following:
if true
% code
function [f,f_time] = Frequency_Calculation_Function(Singal,time)
% This function will locate the zero crossing in your input signal and
% calculate the frequency of that input signal.
% Author: Aubai Al Khatib datum: 19.09.2013
Zeros = Zeros_finding(Singal,time);
f = 1./(2*diff(Zeros));
f_time = Zeros(1:length(f));
end
function [xc] = Zeros_finding(y,x) n = length(y); ind = 1:(n-1);
% list of intervals with a zero crossing
k = find(((y(ind)<=0) & (y(ind+1)>0)) | ...
((y(ind)>=0) & (y(ind+1)<0)));
% list of zero crossings
xc = [];
% intervals where y is zero at both ends of the
% interval are exactly zero are indeterminate.
% the solution may be anywhere on that interval.
% I'll choose to return both endpoints of the
% interval.
L = (y(k)==0) & (y(k+1)==0);
if any(L)
xc = x(k(L));
k(L)=[];
end
% interpolate to find x. I've already removed
% the constant intervals at zero
if ~isempty(k)
s = (y(k+1)-y(k))./(x(k+1)-x(k));
xc = [xc,x(k) - y(k)./s];
end
% patch for last element exactly zero
if y(end)==0
xc(end+1) = x(end);
end
if xc(1) == 0
xc = xc(2:end);
end
end
end
this is the easy version of the code where i show how i calcualted the frequency.
but in the following code i show the complete calculation of the sync. sine wave genertated from that calculated frequency:
if true
% code
function [f,f_time,f_interp,time_Zeros,Wt,Sine,Cose,Signal_Zeros,Signal_Limits,PLL_time] = PLL_final_cal(Signal,time,tol,option)
% This function will locate the zero crossing in your input signal and
% calculate the frequency of that input signal.
% This function will also generate your input signal 'Wt' function and
% regenerate an pu sin & cos signal just like an PLL_Block in simulink.
% *************************************************************************
% the needed inputs:
% 1- Signal: is your input y-axi data (in case of sine wave those are the
% Amplitude)
% 2- time: is your input x-axi data which is the time space of the input
% signal
% 3- tol: is the tolerance in which zeros points will be selected. and it
% should be given in Hz.
% for example if tol = 1 [Hz] is means any Zeros(i+1) point which distance
% from the point before Zeros(i) is in this tolerance band of
% [1/(fn-tol),1/fn,1/(fn+tol)] in S will be accepted as the next frequency
% of f(i+1).
% 4- option: this variable will allow you to choose between doing the PLL
% evaluation according to the batlab simulink block (option = 1) or to
% generate the Wt signal without intgration (option ~= 1)
% both methods should give the same results yet the second one (option ~=1)
% is mor accurate one as Wt is directly built from the linearly
% interpolated zerocorssing function before
% *************************************************************************
% the available outputs:
% 1- f: is the frequency of the input Signal calculated at every Zeros(i)
% points which means f length is the same as the number of Zeros found in
% the input Signal divided by 2 (in the Sine wave the half wave part is
% neglected). mins 1 (due to the "diff" function effect).
% 2- f_time: is the time signals which can be used to "plot" this f signal
% and returns the times where the Zeros has occured and the frequency were
% calculated.
% f_interp & time_Zeros: are the interpulation of the previous signals in
% order to get the output frequency with the same size of the input Signal.
% 3- Wt, Sine & Cose: are the same output of the PLL function block in
% simulink. which means the are the pu represntation of the input Signal
% where there is no harmonics. yet have the same frequency of the input
% Signal.
% 4- Signal_Zeros: is the cut version of the input Signal (to synchronize
% on the first or second Zero crossing)
% 5- Signal_Limits: is the indx representation of the input signals with
% only two values [Start_Value End_Value], in which the complete output was
% created.
% this function file will need the following function files in order to
% work:
% 1- Fast_Harmonics.m: is one is used to extract the fundamental frequency
% out of the input signal.
% 2- Zeros_finding_Rowdata_final.m: this function will be used to extract
% the zero crossing of the input signal and then internally takes only the
% zeros points which fulfill the 'tol' argument mentioned before. so for
% example if you have an input signal where too much harmonics are inside
% and some missleading zero crossing in there is function will try to
% remove those unwanted zeros points.
% Author: Aubai Al Khatib
% 15.01.2014 Created
% 28.01.2013 Modified
Delta_Time_Orginal = time(2)- time(1);
F_Sample = abs(1/Delta_Time_Orginal);
%------------------------------------------------------------------------------------------------------
%---------------------- The Harmonics Spectrum of the complet input signal ----------------------------
%------------------------------------------------------------------------------------------------------
Tp = length(time);
T_S = abs(Tp*Delta_Time_Orginal);
[nyquist_number_I1,Bezugsfrequenz_I1,Signal_amp_I1] = Fast_Harmonics (Signal,F_Sample,T_S);
%plot([0:nyquist_number_I1]*Bezugsfrequenz_I1,Signal_amp_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
F_Spec = [0:nyquist_number_I1]*Bezugsfrequenz_I1; %#ok<*NBRAK>
Fn = round(F_Spec(Signal_amp_I1 == max(Signal_amp_I1)));
[~,Zeros_f,Zeros_f_indx_final] = Zeros_finding_Rowdata_final(Signal,time,Fn,tol);
%[Zeros,k] = Zeros_finding(Signal,time);
Zeros = Zeros_f;
k = Zeros_f_indx_final;
f = 1./(2*diff(Zeros));
f_time = Zeros(1:length(f));
Test_indx = k(1)+2;
if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
signal_test = Signal;time_test = time;
signal_test(k) = 0;time_test(k) = Zeros;
if 0 == 1
Signal_Zeros = Signal(k(2):k(end)-1);
time_Zeros = time(k(2):k(end)-1);
Signal_Final = Signal_Zeros;
Signal_Final(1) = 0; Signal_Final(end) = 0;
Signal_Zeros = Signal_Final;
time_Final = time_Zeros;
time_Final(1) = Zeros(2);time_Final(end) = Zeros(end);
f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f(2:end), diff(k(2:end)),'uniformoutput',0));
time_f_interp = time_Final;
else
Signal_Zeros = signal_test(k(2):k(end)-1);
time_Zeros = time_test(k(2):k(end)-1);
time_Final = time_Zeros;
f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f(2:end), diff(k(2:end)),'uniformoutput',0));
time_f_interp = time_Final;
end
else
signal_test = Signal;time_test = time;
signal_test(k) = 0;time_test(k) = Zeros;
if 0 == 1
Signal_Zeros = Signal(k(1):k(end)-1);
time_Zeros = time(k(1):k(end)-1);
Signal_Final = Signal_Zeros;
Signal_Final(1) = 0; Signal_Final(end) = 0;
Signal_Zeros = Signal_Final;
time_Final = time_Zeros;
time_Final(1) = Zeros(1);time_Final(end) = Zeros(end);
f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
time_f_interp = time_Final;
else
Signal_Zeros = signal_test(k(1):k(end)-1);
time_Zeros = time_test(k(1):k(end)-1);
time_Final = time_Zeros;
f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
time_f_interp = time_Final;
end
end
%f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
%time_f_interp = time_Final;
%time_f_interp = time(k(1):k(end)-1);
time_Zeros_1 = time_Zeros;
time_Zeros = time(k(1):k(end)-1);%time_f_interp;
if option
W = 2*pi*f_interp;
Wt = cumtrapzt(time_Zeros,W);
Wt = mod(Wt,2*pi);
if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
Start_Indx = k(2);
End_Indx = k(end)-1;
Sine = -sin(Wt);
Cose = -cos(Wt);
else
Start_Indx = k(1);
End_Indx = k(end)-1;
Sine = sin(Wt);
Cose = cos(Wt);
end
Signal_Limits = [Start_Indx, End_Indx];
PLL_time = time_f_interp;
else
Zeros_one_period = Zeros(1:2:end);
D_Zeros = kron(Zeros_one_period,[1 1])';
Test = 2:2:length(D_Zeros);
D_Zeros(Test) = D_Zeros(Test)+0.00000001;
A_max = 0;
A_min = 2*pi;
A = [A_min; A_max];%*(pi/180);
R_Z_G = repmat(A,length(D_Zeros),1);
R_Z_G_1 = R_Z_G(1:length(D_Zeros));
%R_Z_G_2 = R_Z_G(1:length(Zeros));
R_Z_G_interp = interp1(D_Zeros,R_Z_G_1,time,'linear');
%R_Z_G_interp = interp1(R_Z_G_1,time,'linear');
R_Z_G_interp_test = R_Z_G_interp(~isnan(R_Z_G_interp));
Wt = R_Z_G_interp_test;
Test_indx = k(1)+2;
if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
Sine = -sin(R_Z_G_interp_test);
Cose = -cos(R_Z_G_interp_test);
else
Sine = sin(R_Z_G_interp_test);
Cose = cos(R_Z_G_interp_test);
end
Time_test = time(isnan(R_Z_G_interp) == 0);
PLL_time = Time_test;
Start_Indx = find(time == Time_test(1));
End_Indx = find(time == Time_test(end));
Signal_Limits = [Start_Indx, End_Indx];
Signal_Zeros = Signal(isnan(R_Z_G_interp) == 0);
end
plots = 0;
if plots
if 0 == 1
figure(); %#ok<*UNRCH>
if 1 == 1
x_min = 1;
x_max = x_min + .01;%time(end);%0.1;
elseif 1 == 0
x_min = time_Zeros(end)-6.5;
x_max = x_min + 0.1;%time(end);%0.1;
elseif 0 == 1
x_min = time_Zeros(1);
x_max = x_min + 0.1;%time(end);%0.1;
else
x_min = time_Zeros(1);
x_max = time_Zeros(end);%time(end);%0.1;
end
subplot(4,1,1);
[AX,ha2,ha3] = plotyy(time,Signal,time_Zeros,Sine);grid on;
set(AX(1),'xlim',[x_min x_max])
set(AX(2),'xlim',[x_min x_max])
set(AX(1),'ylim',[-max(Signal) max(Signal)])
subplot(4,1,2);
plot(time_Zeros,f_interp);grid on;xlim([x_min x_max]);
subplot(4,1,3);
plot(time_Zeros,Wt);grid on;xlim([x_min x_max]);
subplot(4,1,4);
plot(time_Zeros,Sine);grid on;xlim([x_min x_max]);hold on
plot(time_Zeros,Cose,'r');hold off
else
figure(); %#ok<*UNRCH>
if 1 == 1
x_min = 80;
x_max = x_min + 1;%time(end);%0.1;
elseif 1 == 0
x_min = time_Zeros(end)-6.5;
x_max = x_min + 0.5;%time(end);%0.1;
elseif 1 == 0
x_min = time_Zeros(1);
x_max = x_min + 0.1;%time(end);%0.1;
elseif 1 == 0
x_min = time_Zeros(end)-5;
x_max = x_min + 0.1;%time(end);%0.1;
else
x_min = time_Zeros(1);
x_max = time_Zeros(end);%time(end);%0.1;
end
subplot(4,1,1);
[AX,ha2,ha3] = plotyy(time,Signal,time_Final,Sine);grid on;
set(AX(1),'xlim',[x_min x_max])
set(AX(2),'xlim',[x_min x_max])
set(AX(1),'ylim',[-max(Signal) max(Signal)])
subplot(4,1,2);
plot(time_Final,f_interp);grid on;xlim([x_min x_max]);
subplot(4,1,3);
plot(time_Final,Wt);grid on;xlim([x_min x_max]);
subplot(4,1,4);
plot(time_Final,Sine);grid on;xlim([x_min x_max]);hold on
plot(time_Final,Cose,'r');hold off
end
end
end
end

Sign in to comment.

Accepted Answer

Daniel kiracofe
Daniel kiracofe on 13 Aug 2014
Edited: Daniel kiracofe on 13 Aug 2014
the short answer is that you cannot just multiply f by t to get what you want. You need to integrate f with respect to time to get theta. i.e. instead of sin(2 * pi * f * t) you want sin(2 * pi * cumtrapz(t, f)).
this short tutorial I wrote goes into more detail: http://mechanicalvibration.com/Generating_signal_from_freq.html
  2 Comments
Daniel kiracofe
Daniel kiracofe on 14 Aug 2014
I think the easiest way is to just use interp1(). i.e.
Wtnew = interp1( time_Zeros, W, t); Sine = sin(Wtnew);
in this way, Sine will have a sample at every value of the vector t (which could have any sampling rate you want, even a non-constant one). you could also use cubic spline if linear interpolation is not good enough.
note the mod(Wt, 2*pi) is not necessary. it doesn't hurt, but it also is not doing anything for you.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by