# AM.r
# f1: 輸送波の周波数 (Hz)
# f2: 信号波の周波数 (Hz)
# sec: 解析対象時間 (秒)
# n_sample: 解析対象時間の分割数
# DC_offset: AMにおける直流成分
# phase_offset: In-phase信号と輸送波の位相
AM <- function(f1=10000, f2=719, sec=0.1, n_sample=65536, DC_offset=1, phase_offset=pi * 0.3) {
	# プロットエリアを3x3に分割
	pp <- par(mfrow=c(3, 3), mar=c(3,3,3,1))
	on.exit(par(pp))

	# x軸、secで指定された時間を、n_sampleで指定されたサンプル数で分割
	x1 <- seq(0, sec, len=n_sample)
	# 周波数軸
	f <- seq(0, (n_sample - 1) / sec, len=n_sample)

	# 周波数f1およびf2で指定された正弦波、f1が輸送波、f2が信号
	w1 <- sin(x1 * f1 * 2 * pi)
	w2 <- sin(x1 * f2 * 2 * pi) + DC_offset
	# AM信号の作成
	w <- w1 * w2

	# w3がin-phase (I信号)、w4がquadrature (Q信号)生成用の局所発振信号
	w3 <- sin(x1 * f1 * 2 * pi + phase_offset)
	w4 <- sin(x1 * f1 * 2 * pi + phase_offset - pi/2)
	# I/Q抽出に用いる矩形波の作成
	u <- ifelse(w3 > 0, 1, 0)
	v <- ifelse(w4 > 0, 1, 0)

	# 黒: AM波、赤: I信号、青: Q信号
	plot(x1, w, ty="l", xlim=c(0, 0.0005), main="(a)")
	points(x1, u, ty="l", col=2)
	points(x1, v, ty="l", col=4)

	# AM復調用ローパスフィルタの作成、一応信号周波数の2倍まで取る
	lpf <- numeric(length(x1))
	th <- round(f2 * sec * 2)
	lpf[1:th] <- 1
	lpf[1:th + length(x1) - th] <- 1

	# LPF適用前の各信号のFFT。黒: AM波、赤: I信号、青: Q信号
	plot(f, 20*log10(abs(fft(w) / length(w))), ty="l", xlim=c(0, f1 * 2), ylim=c(-80, 0), main="(b)")
	points(f, 20*log10(abs(fft(w * u) / length(w))), ty="l", col=2)
	points(f, 20*log10(abs(fft(w * v) / length(w))), ty="l", col=4)

	# LPF適用後の各信号のFFT。黒: AM波、赤: I信号、青: Q信号
	plot(f, 20*log10(abs(fft(w)) / length(w)), ty="l", xlim=c(0, f1 * 2), ylim=c(-80, 0), main="(c)")
	points(f, 20*log10(abs(fft(w * u) * lpf / length(w))), ty="l", col=2)
	points(f, 20*log10(abs(fft(w * v) * lpf / length(w))), ty="l", col=4)

	# LPF適用後のI信号 (u1)とQ信号 (v1)
	u1 <- Re(fft(lpf * fft(u*w), inv=T)) / length(x1)
	v1 <- Re(fft(lpf * fft(v*w), inv=T)) / length(x1)
	# IとQのプロット
	plot(u1, v1, ty="l", xlim=c(-1,1), ylim=c(-1,1), col=3, main="(d)")
	points(mean(u1), mean(v1))
	abline(h=0, v=0)
	# IとQ、最初の5ms分、I信号が赤、Q信号が青)
	plot(x1, u1, ty="l", xlim=c(0, 0.005), ylim=c(min(c(u1, v1)), max(c(u1, v1))), col=2, main="(e)")
	points(x1, v1, ty="l", col=4)

	# I/Q信号から、輸送波とI信号の位相を推定
	c1 <- mean(u1)
	s1 <- mean(v1)
	c2 <- c1 / sqrt(c1^2 + s1^2)
	s2 <- s1 / sqrt(c1^2 + s1^2)
	# 逆の回転変換を行い、データをX軸上に載せる
	m <- matrix(c(c2, -s2, s2, c2), 2)
	uv <- t(m %*% rbind(u1, v1))
	plot(uv, ty="l", xlim=c(-1,1), ylim=c(-1,1), col=3, main="(f)")
	abline(h=0, v=0)

	# 元の信号と合わせてみる
	# 大体1/3くらいとなる
	plot(x1, w2, ty="l", xlim=c(0, 0.005), ylim=c(0, 4), main="(g)")
	points(x1, uv[,1] * 3, ty="l", col=2)
}

> source("AM.r")
> AM(phase=0)

# 局所発振信号を60度進ませる
> AM(phase=pi/3)