Five Methods for Discrete Spectrum Correction (MATLAB Code)

Five Methods for Discrete Spectrum Correction (MATLAB Code)

1Ratio Correction Method

close all;

clear all;

clc

Fs=1024;

N=1024;

t =(0:N-1)/Fs;

tt=0:N-1;

hanning=0.5-0.5*cos(2*pi*tt/N);

% Hanning window expression

x=5.3*cos(2*pi*253.5453*t+240*pi/180);

% Function to be analyzed

x_hanning=5.3*cos(2*pi*253.5453*t+240*pi/180).*hanning;

% Apply Hanning window

y = fft(x);

yh=fft(x_hanning);

Y=y(1:N/2)/N*2;

Yh=yh(1:N/2)/N*2;

f = Fs/2*linspace(0,1,N/2);

A=abs(Y);

Ah=abs(Yh);

subplot(221);

stem(f,A(1:N/2));

title(‘Uncorrected with Rectangular Window’);

subplot(222);

stem(f,2*Ah(1:N/2));

title(‘Uncorrected with Hanning Window’);

% The factor of 2 in the plot function is the amplitude recovery coefficient

% Apply Rectangular Window

[Amax,k]=max(A);

phmax=angle(Y(k));

if A(k-1)>A(k+1);

deltf=1/(1+A(k)/A(k-1));

f0=(k-1-deltf)*Fs/N

% Corrected frequency

am=Amax/sinc(deltf)

% Corrected amplitude

ph=(phmax+pi*deltf)*180/pi

% Corrected phase

else A(k-1)<A(k+1);

deltf=1/(1+A(k)/A(k+1));

f0=(k-1+deltf)*Fs/N

am=Amax/sinc(deltf)

ph=(phmax-pi*deltf)*180/pi

end

A(k)=am;f(k)=f0;

subplot(223);

stem(f,A(1:N/2));

title(‘Corrected with Rectangular Window’);

% Apply Hanning Window

[Amaxh,kh]=max(Ah);

phmaxh=angle(Yh(kh));

if Ah(kh-1)>Ah(kh+1);

deltfh=(Ah(kh)/Ah(kh-1)-2)/(1+Ah(kh)/Ah(kh-1));

f0h=(kh-1+deltfh)*Fs/N

amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)

phh=(phmaxh-pi*deltfh)*180/pi

else Ah(kh-1)<Ah(kh+1);

deltfh=(Ah(kh)/Ah(kh+1)-2)/(1+Ah(kh)/Ah(kh+1));

f0h=(kh-1-deltfh)*Fs/N

amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)

phh=(phmaxh+pi*deltfh)*180/pi

end

Ah(kh)=amh;

f(kh)=f0h;

subplot(224);

stem(f,Ah(1:N/2));

title(‘Corrected with Hanning Window’);

2Centroid Energy Correction Method

clear all;

clc;

n=input(‘Please enter the number of points for centroid energy correction’)

Fs=1024;N=1024;

t =(0:N-1)/Fs;

tt=0:N-1;

hanning=0.5-0.5*cos(2*pi*tt/N);

x=5.3*cos(2*pi*123.4*t+20*pi/180);

x_hanning=5.3*cos(2*pi*123.4*t+20*pi/180).*hanning;

y = fft(x);

yh = fft(x_hanning);

% The “h” added later indicates Hanning window

Y=y(1:N/2)/N*2;

Yh=yh(1:N/2)/N*2;

f = Fs/2*linspace(0,1,N/2);

A=abs(Y);

Ah=abs(Yh);

subplot(221);

stem(f,A(1:N/2));

title(‘Uncorrected with Rectangular Window’);

grid on;

subplot(222);

stem(f,2*Ah(1:N/2));

title(‘Uncorrected with Hanning Window’);

grid on

% The factor of 2 in the plot function is the amplitude recovery coefficient

G=A.^2;

Gh=Ah.^2;

% Apply Rectangular Window

[Gmax,k]=max(A);

phmax=angle(Y(k));

f0fz=0;

f0fm=0;

for i=-n:n

f0fz=f0fz+(k+i)*G(k+i);

f0fm=f0fm+G(k+i);

end

f0=((f0fz/f0fm)-1)*Fs/N

am=sqrt(1*(f0fm))

% Amplitude recovery coefficient for rectangular window is 1

ph=(phmax+pi*(k-1-f0))*180/pi

A(k)=am;

f(k)=f0;

subplot(223);

stem(f,A(1:N/2));

title(‘Corrected with Rectangular Window’);

grid on

% Apply Hanning Window

[Gmaxh,kh]=max(Ah);

phmaxh=angle(Yh(kh));

f0hfz=0;

f0hfm=0;

for i=-n:n

f0hfz=f0hfz+(kh+i)*Gh(kh+i);

f0hfm=f0hfm+Gh(kh+i);

f0h=f0hfz/f0hfm;

end

f0h=(f0h-1)*Fs/N

amh=sqrt(2.667*(f0hfm))

% 2.66 is the amplitude recovery coefficient for Hanning window

phh=(phmaxh+pi*(kh-1-f0h))*180/pi

Ah(kh)=amh;

f(kh)=f0h;

subplot(224);

stem(f,Ah(1:N/2));

title(‘Corrected with Hanning Window’);

3FFT + FT Spectrum Correction Method

clear all;

clc;

Fs=1024;

N=1024;

t =(0:N-1)/Fs;

tt=0:N-1;

hanning=0.5-0.5*cos(2*pi*tt/N);

windowtype=input(‘Please select window type 1.Rectangular Window 2.Hanning Window’);

if windowtype==1

x=4.2366*cos(2*pi*63.2*t+23.8*pi/180);

% Time domain simulation function

elseif windowtype==2

x=4.2366*cos(2*pi*63.2*t+23.8*pi/180).*hanning;

else

error(‘Invalid selection, please reselect’);

end

y = fft(x);

Y=y(1:N/2+1)/N*2;

f=(0:N/2)*Fs/N;

subplot(211);

plot(f,abs(Y(1:N/2+1)));

grid on

A=abs(Y);

[Amax,k]=max(A);

if windowtype==1

Amax_uncorrect=Amax

% Amplitude before correction

elseif windowtype==2

Amax_uncorrect=Amax*2

end

phmax_uncorrect=angle(Y(k))*180/pi

% Phase before correction

f_uncorrect=k-1

% Frequency before correction

L=80;

% Number of points to refine

deltf=((k+1-1)*Fs/N-(k-1-1)*Fs/N)/L;

% Determine the frequency resolution after refinement

YY=zeros(1,N);

ff=(k-1-1)*Fs/N;

for i=1:L

for ii=1:N

YY(i)=YY(i)+x(ii)*exp(-j*2*pi*ff*(ii-1)/Fs);

end

ff=ff+deltf;

end

if windowtype==1;

A_correct=max(abs(YY)/N*2)

% Amplitude after correction

[YYmax,YYk]=max(YY);

f_correct=(k-1-1)*Fs/N+(YYk-1)*deltf

% Frequency after correction

phmax_correct=angle(YYmax)*180/pi

% Phase after correction

Y(k)=A_correct;f(k)=f_correct;

subplot(212);

plot(f,abs(Y(1:N/2+1)));

grid on

else windowtype==2;

A_correct=2*max(abs(YY)/N*2)

% Amplitude after correction

[YYmax,YYk]=max(YY);

f_correct=(k-1-1)*Fs/N+(YYk-1)*deltf

% Frequency after correction

phmax_correct=angle(YYmax)*180/pi

% Phase after correction

Y(k)=A_correct/2;

f(k)=f_correct;

% Amplitude divided by 2 to ensure the internal values are original values

subplot(212);

plot(f,2*abs(Y(1:N/2+1)));

grid on

% The factor of 2 is the amplitude recovery coefficient

end

~~~~~~~~~~~~~~~~

Running Results:

Theoretical Value

Amplitude: 4.2366 Frequency: 63.2 Phase: 23.8 degrees

With Rectangular Window

Before Correction

Amax_uncorrect =3.962471854470429

phmax_uncorrect =59.672485657326710

f_uncorrect =63

After Correction

A_correct =4.235275753915142

f_correct =63.200000000000003

phmax_correct =23.660397117658885

With Hanning Window

Before Correction

Amax_uncorrect = 4.128431340024557

phmax_uncorrect =59.800004751788116

f_uncorrect =63

After Correction

A_correct =4.236600313219102

f_correct =63.200000000000003

phmax_correct =23.800007456676489

4Time Domain Translation Method

% Fs: Sampling frequency

% N: Number of points for spectrum

% L: Number of translation points

clear all;

clc;

Fs=1024;

N=1024;

L=100;

t =(0:N+L-1)/Fs;

tt=0:N-1;

hanning=0.5-0.5*cos(2*pi*tt/N);

windowtype=input(‘Please select window type 1.Rectangular Window 2.Hanning Window’);

x=10.343*cos(2*pi*298.30453*t+240*pi/180);%.*hanning;

% Time series of L+N points

% Extract two segments of time series

for i=1:N

x1(i)=x(i);

x2(i)=x(i+L);

end

if windowtype==1

y1=fft(x1);

% Perform FFT transformation on the first segment with rectangular window

y2=fft(x2);

% Perform FFT transformation on the second segment with rectangular window

elseif windowtype==2

y1=fft(x1.*hanning);

% Perform FFT transformation on the first segment with Hanning window

y2=fft(x2.*hanning);

% Perform FFT transformation on the second segment with Hanning window

else

error(‘Invalid selection, please reselect’);

end

Y1=abs(y1(1:N/2)/N*2);

% Amplitude of the first segment

Y2=abs(y2(1:N/2)/N*2);

% Amplitude of the second segment

f=(1:N/2)*Fs/N;

subplot(211);

if windowtype==1

plot(f,Y1);

xlabel(‘f’);

ylabel(‘A’);

title(‘Before Correction with Rectangular Window’);

grid on

elseif windowtype==2

plot(f,2*Y1);

xlabel(‘f’);

ylabel(‘A’);

title(‘Before Correction with Hanning Window’);

grid on

end

[Y1Amax,k1]=max(Y1);

[Y2Amax,k2]=max(Y2);

phase1=angle(y1(k1));

% Phase corresponding to the peak of the first segment

phase2=angle(y2(k2));

% Phase corresponding to the peak of the second segment

if windowtype==1

A_uncorrect=Y1Amax

elseif windowtype==2

A_uncorrect=Y1Amax*2

end

f_uncorrect=(k1-1)*Fs/N

% Frequency before correction

phase_uncorrect=phase1*180/pi

% Phase before correction

delt=mod(phase2-phase1,2*pi);

% Adjust delt to be between (-pi, pi)

if delt<-pi

delt1=delt+2*pi;

elseif delt>pi

delt1=delt-2*pi;

else delt1=delt;

end

deltf=delt1/(2*pi*L/N);

f_correct=(k1-1+deltf)*Fs/N

phase_correct=(phase1-deltf*pi)*180/pi

Y_correct=zeros(1,N/2);

if windowtype==1

A_correct=Y1Amax/sinc(deltf)

Y_correct(k1)=A_correct;

f(k1)=f_correct;

subplot(212);

stem(f,Y_correct);

grid on

xlabel(‘f’);

ylabel(‘A’);

title(‘After Correction with Rectangular Window’);

elseif windowtype==2

A_correct=2/sinc(deltf)*Y1Amax*(1-deltf^2)

Y_correct(k2)=A_correct;

f(k2)=f_correct;

subplot(212);

stem(f,Y_correct);

grid on

xlabel(‘f’);

ylabel(‘A’);

title(‘After Correction with Hanning Window’);

end

~~~~~~~~~~~~~~~~

Running Results:

Theoretical Value:

Amplitude: 10.343

Frequency: 298.30453

Phase: 240

With Rectangular Window

Before Correction

A_uncorrect =8.830978191997250

f_uncorrect =298

phase_uncorrect =-65.283647874590685

After Correction

A_correct =10.362781259858098

f_correct =2.983068229594135e+002

phase_correct =-1.205117805690294e+002

With Hanning Window

Before Correction

A_uncorrect =9.739025234720254

f_uncorrect =298

phase_uncorrect =-65.184599901407665

After Correction

A_correct =10.342999982203294

f_correct =2.983045299993504e+002

phase_correct =-1.199999997844760e+002

5Phase Difference Method by Changing Window Length

% Fs: Sampling frequency

% N: Number of points for spectrum

% L: Number of translation points

clear all;

clc;

Fs=1024;

N=1024;

t =(0:N-1)/Fs;

windowtype=input(‘Please select window type 1.Rectangular Window 2.Hanning Window’);

x=10.343*cos(2*pi*298.30453*t+135*pi/180);

%.*hanning;% Time series of L+N points

if windowtype==1

y1=fft(x,N);

% Perform N-point FFT transformation on the signal

y2=fft(x,N/2);

% Perform N/2-point FFT transformation on the signal

elseif windowtype==2

y1=fft(x.*hann(N)’);

% Perform N-point FFT transformation on the signal

y2=fft(x(1:N/2).*hann(N/2)’);

% Perform N/2-point FFT transformation on the signal

else

error(‘Invalid selection, please reselect’);

end

Y1=abs(y1(1:N/2)/N*2);

% Amplitude of the first segment

Y2=abs(y2(1:N/4)/N*4);

% Amplitude of the second segment

f=(1:N/2)*Fs/N;

subplot(211);

if windowtype==1

plot(f,Y1);

xlabel(‘f’);

ylabel(‘A’);

title(‘Before Correction with Rectangular Window’);

grid on

elseif windowtype==2

plot(f,2*Y1);

xlabel(‘f’);

ylabel(‘A’);

title(‘Before Correction with Hanning Window’);

grid on

end

[Y1Amax,k1]=max(Y1);

[Y2Amax,k2]=max(Y2);

phase1=angle(y1(k1));

phase2=angle(y2(k2));

if windowtype==1

A_uncorrect=Y1Amax

% Uncorrected amplitude

elseif windowtype==2

A_uncorrect=Y1Amax*2

end

f_uncorrect=(k1-1)*Fs/N

% Uncorrected frequency

phase_uncorrect=phase1*180/pi

% Uncorrected phase angle

delt=mod(phase1-phase2,2*pi);

% Adjust delt to be between (-pi, pi)

if delt<-pi

delt1=delt+2*pi;

elseif delt>pi

delt1=delt-2*pi;

else delt1=delt;

end

deltf=2*(k2-1)-(k1-1)-2*delt1/pi;

f_correct=(k1-1-deltf)*Fs/N

% Corrected frequency

phase_correct=(phase1+deltf*pi)*180/pi

% Corrected phase

Y_correct=zeros(1,N/2);

if windowtype==1

A_correct=Y1Amax/sinc(deltf)

Y_correct(k1)=A_correct;

f(k1)=f_correct;

subplot(212);

stem(f,Y_correct);

grid on

xlabel(‘f’);

ylabel(‘A’);

title(‘After Correction with Rectangular Window’);

elseif windowtype==2

A_correct=2/sinc(deltf)*Y1Amax*(1-deltf^2)

Y_correct(k2)=A_correct;

f(k2)=f_correct;

subplot(212);

stem(f,Y_correct);

grid on

xlabel(‘f’);

ylabel(‘A’);

title(‘After Correction with Hanning Window’);

end

~~~~~~~~~~~~~~~~

Running Results:

Theoretical Value

Frequency: 298.30453

Amplitude: 10.343

Phase: 135

With Rectangular Window

Before Correction

A_uncorrect =8.836556159610320

f_uncorrect =298

phase_uncorrect =-1.701829534720524e+002

After Correction

A_correct =10.346244554887681

f_correct =2.983047477837129e+002

phase_correct =-2.250375545403775e+002

With Hanning Window

Before Correction

A_uncorrect =9.730664267447002

f_uncorrect =298

phase_uncorrect =-1.702381307852827e+002

After Correction

A_correct =10.334120533701167

f_correct =2.983045300071315e+002

phase_correct =-2.250535320689407e+002

The program is for reference only, and this public account is not responsible for the correctness of the program!

Five Methods for Discrete Spectrum Correction (MATLAB Code)

This article is originally by Zhao Yixu, a member of the Sound and Vibration Forum, and the program references the book “Discrete Spectrum Correction Technology” by Ding Kang and related posts by Yang ZJ, the moderator of the Sound and Vibration Forum. The copyright belongs to the original author, and reprinting requires contacting the original author for authorization and indicating the source (Sound and Vibration Forum: vibunion.com or Sound and Vibration Home WeChat Public Account: vibunion).

Related Readings:
ACauses of Errors in Discrete Spectrum Analysis and Discrete Spectrum Correction Technology
BZoomFFT Has Not Achieved Refinement, Let’s Keep Its Original Name
CBasic Principles of Bandpass Fourier Transform and MATLAB Implementation
DMATLAB Example: Calculation of Spectrum, Power Spectrum, and Power Spectral Density

Related posts

  1. Inventory Optimization Control Based on Teaching-Learning Algorithm with MATLAB Code
  2. MATLAB Signal Processing | Fast DCT Algorithm Using FFT (Complete Code Included)
  3. Analysis of Spring-Mass-Damper System Behavior Using MATLAB
  4. Using Copula for Simulating and Optimizing Market Risk Data VaR Analysis in MATLAB
  5. Qing Reading (27): How to Create a 3D Surface Plot Using MATLAB?
  6. Dynamic Analysis of Power Systems Based on MATLAB | IEEE 9 and IEEE 68 Node Systems
  7. VMware Virtual Machine Performance Management and Resource Online Expansion
  8. Multimodal Targeting Chimeras Enable Integrated Immunotherapy
  9. Is It a Good Thing for Chip Manufacturers to Start Selling Modules?
  10. Installing MSS for MATLAB 2024b Simulink
  11. Using NXP S32G Vehicle Network Processor to Address the Rise of Automotive Gateways and New Vehicle Architectures
  12. Introduction to New Projects in Medical Laboratory – Myositis Antibody Spectrum Testing
  13. Fudan Microelectronics 1H25: Leading Domestic FPGA Manufacturer, New Products like FPAI Worth Looking Forward To

Leave a Comment