周波数が時間的に増加していくプログラムについて

13 views (last 30 days)
創 尾崎
創 尾崎 on 9 Dec 2021
Moved: Atsushi Ueno on 17 Aug 2022
周波数が時間に対して一時間数的に増加するサイン波のプログラムを作りたいと思っています。
式にすると、X = sin(2πft) で、周波数 f は f [Hz] = 12*t [s]という一次関数をイメージしています。
ここで、スペクトログラムに表れた周波数が、実際に想定している周波数と異なってしまうのはなぜか知りたいです。
よろしくお願いいたします。
tend = 5; %計算時間
dt = 0.001; %時間刻み
tnumber = 5001; %時間のサンプリング数
%一時間数的に周波数が増加する波形を生成 ( F(Hz)=12t(s) とする)
F = 0; %最初の周波数
F_k = 12; %周波数の傾き
f = zeros(tnumber,1); %周波数を格納する行列
X = zeros(tnumber,1); %波形を格納する行列
f(1,1) = F; %初速を1行目にセット
T = zeros(tnumber,1); %時間の行列
j = 2;
t = dt;
% t < tend の間の数値積分
while t < tend
f(j,1) = F_k*t; %ローラーの速度V[m/s]が一時間数的に増加
X(j,1) = sin(2*pi*f(j,1)*t); %外力Fの行列
T(j,1) = t;
j = j+1;
t = dt+t;
end
%変位Xのスペクトログラム
figure(1);
spectrogram(X,128,120,128,1000,'yaxis');
xlabel('t[s]');
ylabel('Xの周波数[Hz]');
ylim([0 200]);
box on
%想定している周波数グラフ
figure(2);
plot(T,f,'.','MarkerSize',4);
xlabel('t[s]');
ylabel('f[Hz]');
box on
  2 Comments
Atsushi Ueno
Atsushi Ueno on 9 Dec 2021
Moved: Atsushi Ueno on 17 Aug 2022
chirp関数から必要部のみ抜粋する(ph0=0,f0=0,複素数や曲線を除く)と、下記の様になりました。
じゃなくてとなりますが、何故そうなるのかはわかりません。
type chirp % ← calculateChirp()に着目
function y = chirp(t,varargin) %CHIRP Swept-frequency cosine generator. % Y = CHIRP(T) generates samples of a swept-frequency signal at the time % instances defined in array T. By default, the instantaneous frequency % at time 0 is 0, and the instantaneous frequency one second later is 100 % Hz. % % Y = CHIRP(T,F0,T1,F1) generates samples of a linear swept-frequency % signal at the time instances defined in array T. The instantaneous % frequency at time 0 is F0 Hz. The instantaneous frequency F1 is % achieved at time T1. % % Y = CHIRP(T,F0,T1,F1,method) specifies alternate sweep methods. % Available methods are 'linear', 'quadratic', and 'logarithmic'; the % default is 'linear'. For a logarithmic-sweep, F0 must be greater than % or equal to 1e-6 Hz, which is also the default value. % % Y = CHIRP(T,F0,T1,F1,method,PHI) specifies an initial phase PHI in % degrees. By default, PHI = 0. % % Y = CHIRP(T,FO,T1,F1,'quadratic',PHI,'concave') generates samples of a % quadratic swept-frequency chirp signal whose spectrogram is a parabola % with its concavity in the positive frequency axis. % % Y = CHIRP(T,FO,T1,F1,'quadratic',PHI,'convex') generates samples of a % quadratic swept-frequency signal whose spectrogram is a parabola with % its convexity in the positive frequency axis. % % Y = CHIRP(...,sigtype) specifies sigtype as 'real' or 'complex'; the % default is 'real'. When sigtype is set to 'real', CHIRP generates real % chirp signals. When set to 'complex', CHIRP generates complex chirp % signals. % % EXAMPLE 1: Compute the spectrogram of a linear chirp. % t = 0:0.001:2; % 2 s at 1 kHz sample rate % y = chirp(t,0,1,150); % Start at DC, cross 150 Hz at % % t = 1 s % spectrogram(y,256,250,256,1E3); % Display the spectrogram % % EXAMPLE 2: Compute the spectrogram of a quadratic chirp. % t = -2:0.001:2; % +/-2 s at 1 kHz sample rate % y = chirp(t,100,1,200,'q'); % Start at 100 Hz, cross 200 Hz at % % t = 1 s % spectrogram(y,128,120,128,1E3); % Display the spectrogram % % EXAMPLE 3: Compute the spectrogram of a convex quadratic chirp % t = 0:0.001:1; % 1 s at 1 kHz sample rate % fo = 25; % Start at 25 Hz, % f1 = 100; % go up to 100 Hz % y = chirp(t,fo,1,f1,'q',[],'convex'); % spectrogram(y,256,200,256,1000); % Display the spectrogram. % % EXAMPLE 4: Compute the spectrogram of a concave quadratic chirp % t = 0:0.001:1; % 1 s at 1 kHz sample rate % fo = 100; % Start at 100 Hz, % f1 = 25; % go down to 25 Hz % y = chirp(t,fo,1,f1,'q',[],'concave'); % spectrogram(y,256,200,256,1000); % Display the spectrogram. % % EXAMPLE 5: Compute the spectrogram of a logarithmic chirp % t = 0:0.001:10; % 10 s at 1 kHz sample rate % fo = 10; % Start at 10 Hz, % f1 = 400; % go up to 400 Hz % y = chirp(t,fo,10,f1,'logarithmic'); % spectrogram(y,256,200,256,1000); % Display the spectrogram % % EXAMPLE 6: Compute a complex-valued linear chirp % t = 0:0.001:2; % 2 s at 1 kHz sample rate % fo = -50; % Start at -50 Hz, % t1 = 1; % f1 = 200; % cross 200 Hz at t = 1 s % y = chirp(t,fo,t1,f1,'complex'); % spectrogram(y,256,200,256,1000,'centered'); % Display the spectrogram % % See also GAUSPULS, SAWTOOTH, SINC, SQUARE. % Copyright 1988-2021 The MathWorks, Inc. % References: % [1] Agilent 33220A 20 MHz Function/Arbitrary Waveform Generator % Users Guide, Agilent Technologies, March 2002, pg 298 %#codegen narginchk(1,8); % special case for empty input ([]) if isempty(t) y = zeros(0,1); return end validateattributes(t,{'numeric'},{'real','finite'},'chirp','t',1); t = double(t); % Parse inputs, and substitute for defaults: [f0,t1,f1,method,phi,quadtype,isComplex] = parseOptArgs(varargin{:}); % Initializing output variable y = coder.nullcopy(t); % Computing the signal according to method switch method case 'polynomial' % Polynomial chirp y = cos(2*pi*polyval(polyint(f0),t)); case 'linear' % Linear chirp y = calculateChirp(f0,f1,t1,1,t,phi,isComplex); case 'quadratic' % Compute the quadratic chirp output based on quadtype % Compute the quadratic chirp (upsweep or downsweep) for % complex or concave modes. % For the default 'concave-upsweep' and 'convex-downsweep' modes % call calculateChirp without any changes to the input parameters. % For the forced 'convex-upsweep' and 'concave-downsweep' call % calculateChirp with f0 and f1 swapped and t = fliplr(-t) % For 'convex-upsweep' and 'concave-downsweep' modes if ((f0(1) < f1) && strcmp(quadtype,'convex') || ((f0(1) > f1) && strcmp(quadtype,'concave'))) t = fliplr(-t); f0temp = f1; f1temp = f0(1); else f0temp = f0(1); f1temp = f1; end y = calculateChirp(f0temp,f1temp,t1,2,t,phi,isComplex); case 'logarithmic' % Logarithmic chirp if f0 < 1e-6 coder.internal.error('siglib:chirp:InvalidF0','F0'); else tempVector = repmat(f1/f0,size(t,1),size(t,2)).^(t./t1); instPhi = (t1/log(f1/f0)*f0)*(tempVector-1); if (isComplex) y = exp(1i*2*pi*phi/360).*complex(cos(2*pi*(instPhi)), -cos(2*pi*(instPhi+90/360))); else y = cos(2*pi*(instPhi+phi/360)); end end end end %--------------------------------------------------------------------------- function yvalue = calculateChirp(f0,f1,t1,p,t,phi,isComplex) % General function to compute beta and y for both % linear and quadratic modes. % p is the polynomial order (1 for linear and 2 for quadratic) coder.internal.prefer_const(isComplex); beta = (f1-f0).*(t1.^(-p)); if (isComplex) yvalue = exp(1i*2*pi*phi/360).*complex(cos(2*pi*(beta./(1+p)*(t.^(1+p))+f0*t)),-cos(2*pi*(beta./(1+p)*(t.^(1+p))+f0*t+90/360))); else yvalue = cos(2*pi*(beta./(1+p)*(t.^(1+p))+f0*t+phi/360)); end end %--------------------------------------------------------------------------- function [f0,t1,f1,method,phi,quadtype,isComplex] = parseOptArgs(varargin) % PARSEOPTARGS Parse optional input arguments. % Up to 6 optional input arguments, % f0 - Instantaneous frequency at time 0 % t1 - Reference time % f1 - Instantaneous frequency at time t1 % method - Sweep method % phi - Initial phase % quadtype - Spectrogram shape of quadratic chirp % Flag to keep track of coder target isCodegen = ~coder.target('MATLAB'); % Flag to keep track of output signal type isComplex = false; TFcomprl = false; if ( ~ isempty(varargin)) lastArg = varargin{nargin}; if ( ~ isnumeric(lastArg)) && (strcmpi(lastArg,'complex')) TFrc = strcmpi(lastArg,'complex') || strcmpi(lastArg,'real'); isComplex = true; TFcomprl = true; % checks if the output signal type is coder.const() coder.internal.errorIf((~coder.internal.isConst(lastArg) && ~TFrc),... 'siglib:chirp:CoderConstantExpected'); elseif ( ~ isnumeric(lastArg)) && (strcmpi(lastArg,'real')) TFrc = strcmpi(lastArg,'complex') || strcmpi(lastArg,'real'); isComplex = false; TFcomprl = true; coder.internal.errorIf((~coder.internal.isConst(lastArg) && ~TFrc),... 'siglib:chirp:CoderConstantExpected'); end end % Default value depends upon method f0 = 0; % Flag to check input status of f0 isF0Set = false; % Initialize to default values t1 = 1; f1 = 100; method = 'linear'; phi = 0; quadtype = 'null'; if isCodegen coder.varsize('method',[1,11],[0,1]); coder.varsize('quadtype',[1,7],[0,1]); end % Validation and parsing of each input argument for argIndex = 1:nargin if (TFcomprl) && (argIndex == nargin) break; else input = varargin{argIndex}; if ~isempty(input) switch argIndex case 1 validateattributes(input,{'numeric'},{'vector','finite','real'},'chirp','F0'); % % Deprecated functionality if length(input)>1 % % Y = CHIRP(T,P) specifies a polynomial vector P for the % % instantaneous frequency trajectory of the chirp. if isCodegen coder.internal.error('siglib:chirp:PolynomialNotSupported'); else coder.internal.warning('siglib:chirp:syntaxObsolete'); method = 'polynomial'; end end f0 = double(input); isF0Set = true; case 2 validateattributes(input,{'numeric'},{'scalar','real','positive','finite'},'chirp','T1'); t1 = double(input(1)); case 3 validateattributes(input,{'numeric'},{'scalar','real','finite'},'chirp','F1'); f1 = double(input(1)); case 4 if ischar(input) || isStringScalar(input) % Checking for empty string if ~strcmp(input,'') % % Checking string validity % % Valid input strings optionStrings = {'linear','quadratic','logarithmic'}; method = validatestring(input,optionStrings,'chirp'); end else % Invalid Inputs coder.internal.error('siglib:chirp:InvalidInputType','Method','Chirp','character array',class(input)); end case 5 validateattributes(input,{'numeric'},{'scalar','finite','real'},'chirp','PHI'); phi = double(input(1)); case 6 if ischar(input) || isStringScalar(input) % Checking for empty string if ~strcmp(input,'') % Checking string validity % Valid input strings optionStrings = {'concave','convex'}; quadtype = validatestring(input,optionStrings,'chirp'); end else % Invalid Inputs coder.internal.error('siglib:chirp:InvalidInputType','Shape','Chirp','character array',class(input)); end end end end end % Assigning default value to f0 based on method if ~isF0Set && strcmp(method,'logarithmic') f0 = 1e-6; % Minimum output frequency in Agilent Function generator is 1 microHz end % Quadtype is valid only with quadratic method if ~strcmp(quadtype,'null') && ~strcmp(method,'quadratic') coder.internal.error('siglib:chirp:InvalidMode',quadtype, method); elseif strcmp(method,'quadratic') && strcmp(quadtype,'null') % Determine the shape of the quadratic sweep - concave or convex, % if not defined by user if (f1 > f0) quadtype = 'concave'; % Default for upsweep else quadtype = 'convex'; % Default for downsweep. end end end % [EOF] chirp.m % LocalWords: FO fo Agilent quadtype upsweep downsweep PARSEOPTARGS siglib
f0=0; t1=5; t = 0:0.001:t1; f = 12*t;
X1 = sin(2*pi*f.*t); % ←これだと2倍の周波数が表示される
X2 = sin(2*pi*f(end)*(t.^2)/t1/2); % ←意図した表示になる

Sign in to comment.

Accepted Answer

Toru Ikegami
Toru Ikegami on 10 Dec 2021
スペクトログラムに表示されるような瞬間的な周波数(instantaneous frequency,以下 ) は,信号の位相の時間微分で定義されるからというのがお答えになるでしょうか.
 に対して  ですね.
単純なサイン波はωを時間によらない定数として と書けるので,時刻によらずになります.
ωが時間に依存している場合には事情が変わって となるのでとなります.
特にωが時刻に比例している場合()には となるので,各時刻での瞬間的な角周波数は,設定した周波数の2倍となります.
t = 0:0.001:5;
f = 12*t;
X1 = sin(2*pi*f.*t); %
X2 = sin( pi*f.*t); % <-Instantaneous Frequency が設定の二倍になることを見越して位相を1/2に.
spectrogram(X1,128,120,128,1000,'yaxis');
spectrogram(X2,128,120,128,1000,'yaxis');
  1 Comment
創 尾崎
創 尾崎 on 19 Dec 2021
返信遅くなり申し訳ございません。
スペクトログラムに変換する際に注意が必要だったのですね!
理解できました。ありがとうございました!

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!