用C语言实现下列Matlab程序的功能
程序代码：tic;clc;clear;%close;
filedir=[]; % 设置路径
filename='sp_c2e2.wav'; % 设置文件名
fle=[filedir filename];
[x,fs]=wavread(fle); %读入数据文件
x=x(:,1); %取了所有的数据数据 因为数据是按列存储的，
%x=x(1:58000);
x=x/max(abs(x));
lx=length(x); %求数据长度
length_fft=1024; %fft的长度
lf=length_fft/2+1; %fft之后一半的有效数据
R=xcorr(x); % xcorr(A), when A is a vector, is the auto-correlation sequence. 求自相关序列 R的长度是数据长度的（2*lx-1）行1列
R=R(lx:end); %取从58000到结束的数据
a=find(R==max(R(round(fs/600):fix(fs/60))))-1; % round的作用是最近取整，四舍五入，复数的话把实部和虚部单独取整。fix 无论正负，像0最近取整（向下取整）
f0=fs/a; % a为最大峰值频率的点
n=2*ceil(1.5*a)+1; % length of one frame ceil(X) 是天花板函数。正向最近取整。 %% 取整数towards infinity.--无穷
win=hamming(n)/sum(hamming(n)); %归一化窗函数
[xa,tm]=subframe(x,win,0,fs); %分帧 xa为总的数据， tm为一帧的数据 窗长=帧长，且窗函数没有重叠
%% find peaks and frequency, amplitude, phase
% w=linspace(0,pi,lf);
% e=exp(-1i*(-(n-1)/2:(n-1)/2)'*w);
% fftx=e'*xa;
f=linspace(0,fs/2,lf); %给出513个采样频率等差数列
fftx=fft(xa,length_fft); %xa分帧完一帧数据的长度
fftx=fftx(1:513,:);
m=size(xa,2); % number of frames 获取帧数
DD=cell(3,m ); % cell(M,N) or cell([M,N]) is an M-by-N cell array of empty matrices--方阵.产生3行m列的空方阵
for k=1:m
X=fftx(:,k); %计算频谱，取出第k帧数据
% indmax=peak_selection(abs(X),f0,fs,'h');
indmax=peak_selection(abs(X),fs,'p'); %峰值搜索
Fre=f(indmax);
Am=abs(X(indmax));
Ph=angle(X(indmax));
DD{1,k}=Fre*2*pi; %给单元数组的内容赋值
DD{2,k}=Am;
DD{3,k}=Ph;
end
%% frequency matching
FM=cell(3,m-1);
delta=0.4*f0*2*pi; %频率间隔
for k=1:m-1
a1=DD{1,k};a2=DD{1,k+1}; %a1为第k帧的频率
b1=DD{2,k};b2=DD{2,k+1}; %b1为第k帧的幅度
c1=DD{3,k};c2=DD{3,k+1}; %c1为第k帧的相位
A1=[];A2=[];B1=[];B2=[];C1=[];C2=[]; %设置空矩阵
%上面为准备条件
for i=1:length(a1) %第k帧总频率的个数
e=abs(a1(i)-a2); %计算频差 第k帧的一个频率与第k+1帧的每个频率作差
if min(e)<delta
t=find(e==min(e)); %找到比设定小的频率间隔，把间隔赋值给t，t就是最小间隔
if length(t)~=1
t=t(1);
end
if i==length(a1)
A1=[A1,a1(i)];A2=[A2,a2(t)]; % k 帧中的最后一项
B1=[B1,b1(i)];B2=[B2,b2(t)];
C1=[C1,c1(i)];C2=[C2,c2(t)];
a2(t)=[];b2(t)=[];c2(t)=[]; %删除数组的数据
elseif e(t)<abs(a1(i+1)-a2(t)) % 比相邻的频率差小
A1=[A1,a1(i)];A2=[A2,a2(t)]; %把a1（i）放到A1的后面
B1=[B1,b1(i)];B2=[B2,b2(t)];
C1=[C1,c1(i)];C2=[C2,c2(t)];
a2(t)=[];b2(t)=[];c2(t)=[];
elseif t==1
A1=[A1,a1(i)];A2=[A2,a1(i)];% 比相邻的频率差大且k+1帧中没有m-1频率了，定义death
B1=[B1,b1(i)];B2=[B2,0];
C1=[C1,c1(i)];C2=[C2,c1(i)];
elseif abs(a1(i)-a2(t-1))<delta % 比相邻的频率差大且k+1帧中有 m-1频率
A1=[A1,a1(i)];A2=[A2,a2(t-1)];
B1=[B1,b1(i)];B2=[B2,b2(t-1)];
C1=[C1,c1(i)];C2=[C2,c2(t-1)];
a2(t-1)=[];b2(t)=[];c2(t)=[];
else
A1=[A1,a1(i)];A2=[A2,a1(i)]; % 以上情况都不满足 定义death
B1=[B1,b1(i)];B2=[B2,0];
C1=[C1,c1(i)];C2=[C2,c1(i)];
end
else
A1=[A1,a1(i)];A2=[A2,a1(i)]; % define death 死亡
B1=[B1,b1(i)];B2=[B2,0];
C1=[C1,c1(i)];C2=[C2,c1(i)];
end
end
%上面为第k帧的操作
%下面为第k+1帧的操作
%第m+1帧在第m帧出生
for i=1:length(a2)
t=find(A2<a2(i));
if isempty(t) %判断t是否为空
A1=[a2(i),A1];
A2=[a2(i),A2];
B1=[0,B1];
B2=[b2(i),B2];
C1=[c2(i)-a2(i)*length(A2)/fs,C1];
C2=[c2(i),C2];
else
A1=[A1(1:t(end)),a2(i),A1(t(end)+1:end)];
A2=[A2(1:t(end)),a2(i),A2(t(end)+1:end)];
B1=[B1(1:t(end)),0,B1(t(end)+1:end)];
B2=[B2(1:t(end)),b2(i),B2(t(end)+1:end)];
C1=[C1(1:t(end)),c2(i)-a2(i)*length(A2)/fs,C1(t(end)+1:end)];
C2=[C2(1:t(end)),c2(i),C2(t(end)+1:end)];
end
end
FM{1,k}=[A1;A2]; %A2放在A1的下面
FM{2,k}=[B1;B2];
FM{3,k}=[C1;C2];
end
%% interpolation of alimpulates
IM=cell(2,m); %产生2行m列的空矩阵
for i=1:m-1
e=FM{2,i};
p=size(e,2);
ai=zeros(p,n);
for k=1:p
ai(k,:)=e(1,k)+(e(2,k)-e(1,k))/n*(0:n-1);
end
IM{1,i}=ai;
end
%% interpolation of phase
T=n/fs;
ti=linspace(1/fs,n/fs,n);
for i=1:m-1
e=FM{1,i};t=FM{3,i};
p=size(e,2);c=zeros(p,n);
for k=1:p
q=((t(1,k)+e(1,k)*T-t(2,k))+(e(2,k)-e(1,k))*T/2)/2/pi;
M=round(q); %找到最近的整数
alfa=(t(2,k)-t(1,k)-e(1,k)*T+2*pi*M)*3/(T.^2)-(e(2,k)-e(1,k))/T;
beta=(t(2,k)-t(1,k)-e(1,k)*T+2*pi*M)*(-2)/(T.^3)+(e(2,k)-e(1,k))/(T.^2);
c(k,:)=t(1,k)+e(1,k)*ti+alfa*ti.^2+beta*ti.^3;
% d(k,:)=e(1,k)+2*alfa*ti+3*beta*ti.^2;
% g(k,:)=4*alfa+6*beta*ti;
end
IM{2,i}=c;
end
%% synthesis
s1=zeros(n,m);ns1=zeros(1,m*n);
for i=1:m
for k=0:n-1
s1(k+1,i)=DD{2,i}'*cos(k*(DD{1,i}')/fs+DD{3,i}); %求第i帧的第k+1个数据 /fs规一化频率
end
ns1((i-1)*n+1:n*i)=s1(:,i); %取得第i帧数据
end
soundsc(ns1,fs);
wavwrite(ns1/max(abs(ns1)),fs,'sp_c2e2_a.wav');
s=zeros(n,m);
for i=1:m-1
for k=1:n
s(k,i)=IM{1,i}(:,k)'*cos(IM{2,i}(:,k)); %
end
end
for k=0:n-1
s(k+1,m)=DD{2,m}'*cos(k*(DD{1,m}')/fs+DD{3,m});
end
ns=zeros(1,n*m);
for i=1:m
ns((i-1)*n+1:n*i)=s(:,i);
end
soundsc(ns,fs);
wavwrite(ns/max(abs(ns)),fs,'sp_c2e2_b.wav'); %正弦插值合成 可能是最终需要的合成
soundsc(x,fs);
wavwrite(x/max(abs(x)),fs,'sp_c2e2_c.wav'); %直接合成
toc;
%%
% [xa,tm]=subframe(x,rectwin(n),0,fs);
% w=linspace(0,pi,lf);
% f=linspace(0,fs/2,lf);
% e=exp(-1i*(-(n-1)/2:(n-1)/2)'*w);
% fftx=e'*xa;
% ffts=e'*s;
% figure;plot(unwrap(angle(fftx(:,1))))
% figure;plot(unwrap(angle(ffts(:,1))))
