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

周波数が時間に対して一時間数的に増加するサイン波のプログラムを作りたいと思っています。

式にすると、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

[Atsushi Ueno]

スイープ周波数の余弦信号 - MATLAB chirp - MathWorks 日本

>スペクトログラムに表れた周波数が、実際に想定している周波数と異なってしまうのはなぜか知りたいです

理屈が分からず説明出来ませんが、chirp関数を使って同じ目的の信号を作ると意図した通りに表示されました。

t = 0:0.001:5; % f = 12*t; % X = sin(2*pi*f.*t); % ←これだと2倍の周波数が表示される
X = chirp(t,0,5,60); % 0[s]で0[Hz], 5[s]で60[Hz]
spectrogram(X,128,120,128,1000,'yaxis');
xlabel('t[s]'); ylabel('Xの周波数[Hz]'); ylim([0 200]);

[Atsushi Ueno]

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); % ←意図した表示になる

Answer 1

[Toru Ikegami]

スペクトログラムに表示されるような瞬間的な周波数（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');

[創 尾崎]

返信遅くなり申し訳ございません。

スペクトログラムに変換する際に注意が必要だったのですね！

理解できました。ありがとうございました！