|
|
发表于 2015-2-26 23:20:31
|
显示全部楼层
这个问题其实可以表述为：信号为频率已知但幅度和相位未知的正弦波，淹没在干扰和噪声中，现测量到一组值[y(0), .... y(N-1)]，要求根据测量值估计信号的幅度和相位，这个问题可以通过最小二乘来求得。假设信号x(t) = a*cos(2*pi*f*t)+b*sin(2*pi*f*t)，噪声为n(t)，那么测量值为y(t) = x(t) + n(t)，以采样频率fs采样后得到N个采样值y(k) = x(k) + n(k) = a*cos(2* pi * k * f/fs )+b*sin(2*pi*k*f/fs) + n(k)，最小二乘解是
[a b]^T = (A^T * A)^-1 * A^T*Y，其中A = [cos(2* pi * 0 * f/fs ) sin(2* pi * 0 * f/fs );cos(2* pi * 1 * f/fs ) sin(2* pi * 1 * f/fs );......;cos(2* pi * (n-1) * f/fs ) sin(2* pi * (n-1) * f/fs );]，Y=[y(0) y(1) ... y(N-1)]，估计出a,b后，就可以计算出幅度和相位了。
例子：信号频率为45.678Hz，幅度为0.9，相位为10度，干扰信号为50Hz，幅度为10的正弦波，不加噪声。Matlab代码仿真最终估计出来的幅度为0.89846，相位为9.7638度
f = 45.678; % 信号频率
fi = 50; % 干扰信号频率
fs = 150; % 采样率
N = 2048; % 采样点数，采样越多，估计结果越精确
t = [0:N-1]' / fs;
x = 0.9 * sin(2 * pi * f * t + 10 * (pi / 180)); % 原信号
xi = 10 * sin(2 * pi * fi * t); % 干扰信号
xn = 0 * randn(N,1); % 噪声
y = x + (xi + xn); % 被干扰和噪声污染的测量信号
% 最小二乘
A = [cos(2*pi*f*t) sin(2*pi*f*t)]; % 构造矩阵A
result = (A' * A)^-1 * A' * y; % 最小二乘，求出a和b
xp = result(1) * cos(2*pi*f*t) + result(2) * sin(2*pi*f*t);% 根据估计出来的参数重构信号
plot([y, x, xp], '.-'); % 将测量信号，原信号，重构后的信号画在同一幅图中方便比较
amp = sqrt(result(1)^2 + result(2)^2);
phase = atan(result(1) / result(2)) * 180 / pi; |
|
楼主
哥不写
