DFT计算相位差（MATLAB程序）

%本程序通过过给定不同相位差的信号，然后通过DFT分析，计算出相位差，并绘制成输入输出图形 %本程序暂时只研究单一输入频率 %采样频率为输入频率的128倍频。 %本程序主要研究输入信号中含有高次谐波的情况，通过仿真分析，高次谐波几乎不影响基波相位的测量 clc; clear; clf; f0=input('the frequence of input signal:'); fs=f0*128; N=128; %sample 128 dots，representing one complete cycle of input signal. n=0:N-1; %discrete sequence t=n/fs; %discrete the time sequence, so as to plot the time-domain wave xlabeled by "t" f=n*fs; %discrete the physical frequence, so as to plot the frequence response wave xlabeled by "f" %phase difference sequence with 2*pi split into 30 elements; k=12; %k is the theoretical_diff=-pi:2*pi/k:pi; calculated_diff=zeros(1,k+1); %initialize the sequence with 0,modify the elements in the following loops error=zeros(1,k+1); %initialize the sequence with 0,modify the elements in the following loops x1=sin(2*pi*f0*t)+0.2*sin(2*pi*f0*3*t)+0.1*sin(2*pi*f0*5*t); %signal one pulluted by harmonic wave y1=fft(x1,N); %fft the first sequence figure(1); plot(t,x1,'b'); hold on; x2=sin(2*pi*f0*t+pi/6)+0.2*sin(2*pi*f0*7*t); plot(t,x2,'g'); hold off; title('input signals'); xlabel('time t'); ylabel('amplitude y'); legend('signal 1','signal 2'); grid on; ang1=angle(y1); %calculate the phase used later for i=1:(k+1) x2=sin(2*pi*f0*t+theoretical_diff(i))+0.3*sin(2*pi*f0*3*t); %signal two,with a phase difference y2=fft(x2,N); %fft the second sequence ang2=angle(y2); %calculate the phase response ang_diff=ang2-ang1; calculated_diff(i)=ang_diff(2); %ang_diff(2) represents the base frequence,the nomalized w is 2pi/N, the physical frequence is f0. %N=256,512...k*fs,that means k cycles are sampled. then using ang_diff(2k-1) which represnets base frequence. if(calculated_diff(i)>0)&&(theoretical_diff(i)<0) calculated_diff(i)=calculated_diff(i)-2*pi; end if(calculated_diff(i)<0)&&(theoretical_diff(i)>0) calculated_diff(i)=calculated_diff(i)+2*pi; end error(i)=(calculated_diff(i)-theoretical_diff(i))*57.3; i=i; end; figure(2); plot(theoretical_diff/pi,calculated_diff/pi,'*r-'); hold on; plot(theoretical_diff/pi,error,'ob-') hold off; title('phase difference diagram'); xlabel('theoretical phase difference \omega'); ylabel('calculated phase difference,unit（\pi）/error,unit（\circ） ') axis([-1 1 -1 1]); grid on; legend('calculated value','error'); figure(3); stem(n*fs/N,abs(y1)/N);; title('frequence-amplitude diagram'); xlabel('frequence'); ylabel('normalized amplitude'); grid on; axis([0 128*f0 0 1]);