双谱线插值与三谱线插值FFT的MATLAB实现

双谱线插值
f0=50.1; % 基波频率
fs=1500; % 采样频率
N=512; % 数据长度
n=0:N-1; % 数据索引
rad=pi/180; % 角度和弧度的转换因子
xb=[1,0.02,0.1,0.01,0.05,0,0.02,0,0.01]; % 谐波幅值
Q=[-23,115.6,59.3,52.4,123.8,0,-31.8,0,-63.7]rad; % 谐波初始相位
s=zeros(1,N); % 初始化
M=9; % 谐波个数
for i=1:M % 产生谐波信号
s=s+xb(i)cos(2pif0in./fs+Q(i));
end

w=0.5-0.5cos(2pi*n./N); % 海宁窗
x=s.*w; % 信号乘以窗函数
v=fft(x,N); % FFT
u=abs(v); % 取频谱的幅值
k1=zeros(1,M); % 初始化
k2=zeros(1,M);
A=zeros(1,M);
ff=zeros(1,M);
Ph=zeros(1,M);
df=fs/N; % 频率分瓣率

for i=1:M % 计算基波和各阶谐波的参数
if i==1 % 若计算基波,在35-65Hz区间中寻找最大峰值
n1=fix(35/df); n2=fix(65/df); % 求出35Hz和65Hz对应的索引号
else % 若计算谐波,从该谐波理论值-15和+15的区间中寻找最大值
n1=fix((iff(1)-15)/df); % 求出区间对应的索引号
n2=fix((iff(1)+15)/df);
end
[um,ul]=max(u(n1:n2)); % 在区间中找出最大值
k1(i)=ul+n1-1; % 给出最大值的索引号
% 判断峰值在最大值左边还是右边,如果峰值在最大值左边,把k1(i)进行修正
if u(k1(i)-1)>u(k1(i)+1),
k1(i)=k1(i)-1;
end
k2(i)=k1(i)+1; % 求出k2(i),使峰值永远在k1(i)和k2(i)之间
y1=u(k1(i)); % 求出y1和y2
y2=u(k2(i));
b=(y2-y1)/(y2+y1); % 按式(7-2-5)计算出beta
% 海宁窗的beta对alpha的表示式
a=1.5b;
ff(i)=(k1(i)-1+a+0.5)fs/N; % 求出谐波的频率
% 按式(7-2-9)依海宁窗的窗函数关系计算出谐波的幅值
A(i)=(y1+y2)(2.35619403+1.15543682a^{2+0.32607873*a}4+…
0.07891461a^6)/N;
if y1>=y2
Ph(i)=phase(v(k1(i)))+pi/2-pi(a+0.5); % 求出谐波的初始相位
Ph(i)=Ph(i)-(Ph(i)>pi)2pi+(Ph(i)<-pi)2pi; % 对相位进行修正
end
if y1<y2
Ph(i)=phase(v(k1(i)))+pi/2-pi*(a-0.5); % 求出谐波的初始相位
Ph(i)=Ph(i)-(Ph(i)>pi)2pi+(Ph(i)<-pi)2pi; % 对相位进行修正
end
% 若幅值过小设为0,并对频率相位修正
if A(i)<0.0005, A(i)=0; ff(i)=i*ff(1); Ph(i)=0;
end
% 显示谐波参数
fprintf(’%4d %5.6f %5.6f %5.6f\n’,i,A(i),ff(i),Ph(i)/rad);
end

三谱线插值程序：
clc
f0=50.1; % 基波频率
fs=1500; % 采样频率
N=512; % 数据长度
n=0:N-1; % 数据索引
rad=pi/180; % 角度和弧度的转换因子
xb=[1,0.02,0.1,0.01,0.05,0,0.02,0,0.01]; % 谐波幅值
Q=[-23,115.6,59.3,52.4,123.8,0,-31.8,0,-63.7]rad; % 谐波初始相位
s=zeros(1,N); % 初始化
M=9; % 谐波个数
for i=1:M % 产生谐波信号
s=s+xb(i)cos(2pif0in./fs+Q(i));
end

for i=1:M % 计算基波和各阶谐波的参数
if i==1 % 若计算基波,在35-65Hz区间中寻找最大峰值
n1=fix(35/df); n2=fix(65/df); % 求出35Hz和65Hz对应的索引号
else % 若计算谐波,从该谐波理论值-15和+15的区间中寻找最大值
n1=fix((iff(1)-15)/df); % 求出区间对应的索引号
n2=fix((iff(1)+15)/df);
end
[um,ul]=max(u(n1:n2)); % 在区间中找出最大值
k2(i)=ul+n1-1; % 给出最大值的索引号
k1(i)=k2(i)-1;
k3(i)=k2(i)+1;
y1=u(k1(i)); % 求出y1和y2,y3
y2=u(k2(i));
y3=u(k3(i));
b=(y3-y1)/y2; % 按式(7-2-5)计算出beta
% 海宁窗的beta对alpha的表示式
a=0.66666287b-0.07398320b^{3+0.01587358*b}5-0.00311639b^7 ;
ff(i)=(k2(i)-1+a)fs/N; % 求出谐波的频率
% 按式(7-2-9)依海宁窗的窗函数关系计算出谐波的幅值
A(i)=(y1+2y2+y3)(1.33333327+0.52658791a^{2+0.11699742*a}4+…
0.02103885a^6)/N;
Ph(i)=phase(v(k2(i)))+pi/2-api; % 求出谐波的初始相位
Ph(i)=Ph(i)-(Ph(i)>pi)2pi+(Ph(i)<-pi)2pi; % 对相位进行修正
% 若幅值过小设为0,并对频率相位修正
if A(i)<0.0005, A(i)=0; ff(i)=iff(1); Ph(i)=0;
end
% 显示谐波参数
fprintf(’%4d %5.6f %5.6f %5.6f\n’,i,A(i),ff(i),Ph(i)/rad)
end

双谱线插值与三谱线插值FFT的MATLAB实现相关推荐

图形图像处理－之－高质量的快速的图像缩放中篇二次线性插值和三次卷积插值
from:http://blog.csdn.net/housisong/article/details/1452249 图形图像处理-之-高质量的快速的图像缩放中篇二次线性插值和三次卷积插值 ...
MATLAB一维插值和二维插值
插值问题描述:已知一个函数上的若干点,但函数具体表达式未知,现在要利用已知的若干点求在其他点处的函数值,这个过程就是插值的过程. 1.一维插值一维插值就是给出y=f(x)上的点(x1,y1),(x2 ...
python数学建模（三）插值常用库和模块
文章目录 (一)本文数据资料下载 (二)简单介绍一下定义 (三)介绍我们可能用到的模块和代码(重点) 3.1 scipy.interpolate 模块 3.1.1 一维插值函数 (interp1d) ...
iDesktop点数据集构建DEM时三种插值方式的选择
转自:https://blog.csdn.net/supermapsupport/article/details/76252498 点数据构建DEM的时候,可以选择三种插值方式,分别是不规则三角网法, ...
数学建模准备插值（拉格朗日多项式插值，牛顿多项式插值，分段线性插值，分段三次样条插值，分段三次Hermite插值）
文章目录摘要(必看) 0 基础概念什么是插值插值用途什么是拟合插值和拟合的相同点插值和拟合的不同点 1 常用的基本插值方法 1.1 多项式插值法 1.1.1 拉格朗日多项式插值法多项式插 ...
基于神经网络多项式插值的图像超分辨重构研究-附Matlab代码
⭕⭕ 目录 ⭕⭕ ✳️ 一.引言 ✳️ 二.基于单帧图像的超分辨率重构技术 ✳️ 2.1 最近邻域插值法 ✳️ 2.2 双线性插值法 ✳️ 2.3 双三次插值法(Keys'插值) ✳️ 三.神经网络 ...
图像处理中的“内插”是什么？插值、图像内插值、图像间插值、重取样（用已知数据来估计未知位置的数值的处理）（最近邻内插法、双线性内插）
图像插值是在基于模型框架下,从低分辨率图像生成高分辨率图像的过程,用以恢复图像中所丢失的信息.图象插值方法有:最近邻插值,双线性插值,双平方插值,双立方插值以及其他高阶方法. 在很多情况下,人们需要对 ...
插值与拟合 (一) ：拉格朗日多项式插值、Newton插值、分段线性插值、Hermite插值、样条插值、 B 样条函数插值、二维插值
插值:求过已知有限个数据点的近似函数. 拟合:已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义下它在这些点上的总偏差最小. 插值和拟合都是要根据一组数据构造一个函数作为近似,由于近似 ...
一元函数的插值c语言,一元函数插值-Read.doc
一元函数插值目标: 对比插值和数据拟合的思想一元函数插值(线性插值和抛物线插值)及程序实现基本理论插值:采集到的数据可以看作是精确数据,要求函数或曲线通过这些数据点. 常用插值方法:一元函数插 ...

双谱线插值与三谱线插值FFT的MATLAB实现

双谱线插值与三谱线插值FFT的MATLAB实现相关推荐

最新文章

热门文章