20周期ごとに異なる​周波数の正弦波を乗算​するプログラム

現在下記のプログラムを作成しており、繰り返し波形に対し周波数を変化させた正弦波を乗算するプログラムを作りたいのですがのこぎりのような波形になってしまい上手く振幅変調ができません。

%% Read Audio File

[y,Fs]=audioread('piano_ra_2.wav');

dt=1/Fs;

duration=(length(y))*(1/Fs);

%% Cut Waveform

f0=439.8;

period=1/f0;

n0=floor(length(y)/2);

n=round(period/dt);

y_one_period=y([n0:n0+n]);

n_period_1=860;

f_lo=5;

T=0.8*(1/f_lo);

y2=repmat(y_one_period,[n_period_1,1]);

t2=1/Fs:1/Fs:43860*(1/Fs);

y_lo=((1/9)*cos(2*pi*(1/T)*t2))+1;

figure;

plot(t2,y_lo);

set(gca,'FontSize',15);

xlabel('Time (s)','FontSize',25)

ylabel('Relative Amplitude','FontSize',25)

%% rand関数を用いた正規分布の作成（平均値と分散）

va=2;

r=va.*randn(20,1)+6.25;

histogram(r);

%% 20周期単位に周波数変更

r1=6.3849;

r2=5.8758;

r3=6.8335;

r4=8.2254;

r5=7.0359;

r6=6.6391;

r7=6.8096;

r8=6.3524;

r9=4.7011;

r10=7.8236;

r11=9.0678;

r12=5.1818;

r13=10.1055;

r14=5.8975;

r15=5.7625;

r16=4.4548;

r17=4.6653;

r18=4.3441;

r19=6.9578;

r20=9.4441;

a=1/9;

n_period_2=43;

y3=repmat(y_one_period,[n_period_2,1]);

t3=1/Fs:1/Fs:(43860/20)*(1/Fs);

y_r1=a*cos(2*pi*r1*t3)+1;

y_r2=a*cos(2*pi*r2*t3)+1;

y_r3=a*cos(2*pi*r3*t3)+1;

y_r4=a*cos(2*pi*r4*t3)+1;

y_r5=a*cos(2*pi*r5*t3)+1;

y_r6=a*cos(2*pi*r6*t3)+1;

y_r7=a*cos(2*pi*r7*t3)+1;

y_r8=a*cos(2*pi*r8*t3)+1;

y_r9=a*cos(2*pi*r9*t3)+1;

y_r10=a*cos(2*pi*r10*t3)+1;

y_r11=a*cos(2*pi*r11*t3)+1;

y_r12=a*cos(2*pi*r12*t3)+1;

y_r13=a*cos(2*pi*r13*t3)+1;

y_r14=a*cos(2*pi*r14*t3)+1;

y_r15=a*cos(2*pi*r15*t3)+1;

y_r16=a*cos(2*pi*r16*t3)+1;

y_r17=a*cos(2*pi*r17*t3)+1;

y_r18=a*cos(2*pi*r18*t3)+1;

y_r19=a*cos(2*pi*r19*t3)+1;

y_r20=a*cos(2*pi*r20*t3)+1;

figure;

plot(t3,y_r1);

set(gca,'FontSize',15);

xlabel('Time (s)','FontSize',25)

ylabel('Relative Amplitude','FontSize',25)

x_r=[y_r1,y_r2,y_r3,y_r4,y_r5,y_r6,y_r7,y_r8,y_r9,y_r10,y_r11,y_r12,y_r13,y_r14,y_r15,y_r16,y_r17,y_r18,y_r19,y_r20];

y_len=20*length(t3);

y_r=0:y_len-1;

figure;

plot(y_r,x_r);

set(gca,'FontSize',15);

xlabel('Time (s)','FontSize',25)

ylabel('Relative Amplitude','FontSize',25)

t4=1/Fs:1/Fs:43860*(1/Fs);

figure;

plot(t4,y_r);

set(gca,'FontSize',15);

xlabel('Time (s)','FontSize',25)

ylabel('Relative Amplitude','FontSize',25)

%% 混合器

y2=y2.';

mix_r=y2.*y_r;

figure;

plot(t4,mix_r);

set(gca,'FontSize',15);

xlabel('Time (s)','FontSize',25)

ylabel('Relative Amplitude','FontSize',25)

Answer 1

[Naoya]

以下の様なイメージとなりますでしょうか？

(以下は サンプリング周波数 1KHz で 10秒間の信号を作成しています。

1秒毎に ランダムな周波数を与えて、周波数の変わり目の位相も連続的にしている例となります)

Fs = 1000; % サンプリング周波数(Hz)
t = 0:1/Fs:10-1/Fs; % 周波数

cnt = ones(1,10000); % カウンター用ベクトル

% 1秒ごとの周波数を設定
f = [12.4 13.4 7.2 3.3 7.3 1.4 5.4 8.9 9.0 2.2];
f = repmat(f,[1000, 1]);
f = f(:)';

% 正弦波の位相 ベクトル作成
ft = cumsum(cnt.*f)/Fs;

% 周波数を周期的に変化する正弦波作成
sig = sin(2*pi* ft );
plot(t, sig);