Improving personal Fourier Interpolation Model

Hi i'm trying to interpolate a cyclical set of data through a fourier transform model. I can't just use intrpft because i need the coefficients (amplitude, frequency and phase) to project the wave equation into future. I attach the program.

Basically i use the FFT to find n cycles amplitude and frequency. Then i sync each equation with the data. At the end i superimpose all the n cycles to find the wave equation.

I reached an R2 of 0.71.

The model interpolates very well frequencies and phases but there is some problems with amplitudes i guess.

I wonder if somebody can improve precision further.

Thanks

Answer 1

[Star Strider]

It is difficult to understand what you are doing, or the question you aree asking.

I ran your code and it appears to work.

The phase information is straightforward to calculate:

FFT = fft(yv,nfft)/Datapoints;

phase = angle(FFT); % <— ADDED Calculate Phase (radians)

FFT = 2*abs(FFT(1:nfft/2+1));

f = fs/2*linspace(0,1,nfft/2+1);

The ‘coefficients’ in a Fourier transform are the magnitudes of the cosine (real) and sine (imaginary) terms at each frequency.

You can use the sine and cosine terms at each frequency to reconstruct the timne-domain signal.

.

[Giacomo]

I tried to implement your suggestion calculating the phase in this way. However things get worse.

I attach the script.

I'm simply asking if there is a way to improve R2, so that interpolation follows better the data.

I think that the problem is in the amplitude calculation. Frequencies and phases are ok.

I calculate the phase finding the maximum of the correlation function between every single cycle and the data.

[Star Strider]

My approach to reconstructing the time-domaiun waveform from the Fourier coefficients (R² = 0.996) —

x =[0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19    20 21    22    23    24    25    26    27    28    29    30    31    32    33    34    35    36    37    38    39    40    41 42    43    44    45    46    47    48    49    50    51    52    53    54    55    56    57    58    59    60    61    62 63    64    65    66    67    68    69    70    71    72    73    74    75    76    77    78    79    80    81    82    83 84    85    86    87    88    89    90    91    92    93    94    95    96    97    98    99   100   101   102   103   104 105   106   107   108   109   110   111   112   113   114   115   116   117   118   119   120   121   122   123   124   125 126   127   128   129   130   131   132   133   134   135   136   137   138   139   140   141   142   143   144   145   146 147   148   149   150   151   152   153   154   155   156   157   158   159   160   161   162   163   164   165   166   167 168   169   170   171   172   173   174   175   176   177   178   179   180   181   182   183   184   185   186   187   188 189   190   191   192   193   194   195   196   197   198   199   200   201   202   203   204   205   206   207   208   209 210   211   212   213   214   215   216   217   218   219   220   221   222   223   224   225   226   227   228   229   230 231   232   233   234   235   236   237   238   239   240   241   242   243   244   245   246   247   248   249   250   251 252   253   254   255   256   257   258   259   260   261   262   263   264   265   266   267   268   269   270   271   272 273   274   275   276   277   278   279   280   281   282   283   284   285   286   287   288   289   290   291   292   293 294   295   296   297   298];
y =[-10.5109   -1.2195    7.8469   14.1429   16.2819   14.6567   11.2349    8.6033    8.7078   11.8968   16.7308   20.6538   21.1976   17.1386    9.0625   -0.9103  -10.0699  -16.2555  -18.6490  -17.8892  -15.5312  -13.1869  -11.8010  -11.3909  -11.2983  -10.7266   -9.2172   -6.8066   -3.8283   -0.5531    3.0506    7.2480   12.1846   17.4779   22.0492   24.3720   23.0792   17.6471    8.7990   -1.6092  -11.2030  -17.9301  -20.7384  -19.7912  -16.1598  -11.2007   -5.9639   -0.9310    3.8327    8.2740   11.9973   14.2207  14.1022   11.2560    6.1434    0.0715   -5.2581   -8.4492   -9.0052   -7.5011   -5.2318   -3.5181   -3.0353   -3.5248   -4.0378 -3.5722   -1.7391    0.9238    3.1897    3.8204    2.3237   -0.6352   -3.4286   -4.2130   -1.8732    3.3114    9.6375   14.6372    16.1367   13.2253    6.6788   -1.3565   -8.2684  -12.0471  -12.0361   -9.0487   -4.8460   -1.2646    0.5771    0.6378   -0.2435    -0.9455   -0.7268    0.3869    1.6692    2.1536    1.2313   -0.9304   -3.3716   -4.7455   -3.9615   -0.7664    4.0153    8.7097   11.4121   10.7500    6.4734   -0.3733   -7.7794  -13.4745  -15.7392  -14.0051   -9.0118   -2.4672    3.6402    7.7917    9.4438   9.0905    7.8818    6.9889    7.0208    7.7607    8.3378    7.7356    5.3731    1.4572   -3.0786   -7.0221   -9.4841  -10.3331 -10.2179  -10.1477  -10.8409  -12.1942  -13.1898  -12.3401   -8.4894   -1.5794    7.0370   14.9522   19.6282   19.4898   14.6770    7.0773   -0.4123   -5.0824   -5.5302   -2.1934    2.9356    7.2711    8.8090    6.9740    2.7630   -1.8432   -4.8836   -5.3432   -3.5700   -1.0217    0.4950   -0.3003   -3.5852   -8.3959  -13.0926  -16.0629  -16.3249  -13.7819   -9.0856   -3.2546    2.7213    8.1070   12.4336   15.3871   16.7195   16.2816   14.1588   10.7981    6.9895    3.6470    1.4554    0.5552    0.4452    0.1855  -1.1633   -4.0650   -8.2230  -12.6336  -15.9623  -17.0760  -15.4899  -11.5380   -6.2099   -0.7514    3.7744    6.7799    8.2032  8.3496    7.6510    6.4775    5.0813    3.6523    2.4035    1.5944    1.4562    2.0573    3.2000    4.4279    5.1684    4.9524 3.6069    1.3246   -1.4270   -4.1073   -6.3144   -7.8987   -8.9312   -9.5519   -9.7939   -9.4848   -8.2921   -5.8971   -2.2196  2.4193    7.2823   11.3738   13.7237   13.6958   11.2107    6.8034    1.4920   -3.5023   -7.0973   -8.6461   -8.0951   -5.9509 -3.0666   -0.3190    1.7049    2.8917    3.5746    4.2977    5.4674    7.0454    8.4417    8.6881    6.8505    2.5180   -3.8545 -10.9311  -16.8617  -19.9141  -19.1236  -14.6870   -7.9182   -0.7697    4.8922    7.9559    8.4032    7.1844    5.7055    5.1666   6.0469    7.9541    9.8788   10.7135    9.7826    7.1424    3.5299    0.0083   -2.5093   -3.6337   -3.6112   -3.1454   -3.0333 -3.7945   -5.4569   -7.5774   -9.4631  -10.4763  -10.2829   -8.9475   -6.8554   -4.5183   -2.3599   -0.5720    0.9113    2.3268 3.9330    5.8640    8.0423   10.1782   11.8523   12.6467   12.2813   10.7082    8.1341    4.9623    1.6760   -1.2969   -3.6892];
% yv=y';
% Length=x(end);
% Datapoints=length(x);
% fs=round(Datapoints/Length);

fs = 1/(x(2) - x(1))

fs = 1

L = numel(x);
nfft = 2^nextpow2(L)*2;

FFT = fft(y(:),nfft)/L;
% FFT = 2*abs(FFT(1:nfft/2+1));
f = fs/2*linspace(0,1,nfft/2+1);

A = real(FFT(1:numel(f)));
ph = angle(FFT);
for k = 1:numel(f)
    Re(k,:) = A(k)*cos(2*pi*f(k)*x + angle(k));
end

Res = sum(Re)

Res = 1×299

-17.8169 -0.8297 6.9001 12.3234 14.1220 12.7633 9.8008 7.5804 7.6372 10.4003 14.5064 17.8979 18.3308 14.8882 7.9409 -0.5650 -8.4401 -13.7034 -15.7854 -15.1021 -13.1159 -11.0761 -9.9222 -9.5383 -9.4918 -8.9696 -7.7100 -5.6133 -3.0961 -0.2592

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Mdl = fitlm(y, Res)

Mdl = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate       SE        tStat      pValue  
                   ________    _________    ______    _________

    (Intercept)    0.16637      0.030176    5.5135    7.649e-08
    x1             0.85999     0.0032431    265.18            0


Number of observations: 299, Error degrees of freedom: 297
Root Mean Squared Error: 0.521
R-squared: 0.996,  Adjusted R-Squared: 0.996
F-statistic vs. constant model: 7.03e+04, p-value = 0

figure
plot(x, Res, 'LineWidth',1.5, 'DisplayName','Reconstructed Signal')
hold on
plot(x, y, '--', 'LineWidth',1, 'DisplayName','Original Signal')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

[K,X] = ndgrid(1:numel(f),x);

figure
plot3(X, K, Re.')
grid on
xlabel('Time')
ylabel('Term')
title('Re')

return

[y_peak,x_peak]=findpeaks(FFT,f);

%wave equation yn=sum(An*sin(2*pi*Fn*(x+Pn)) n=number of peaks (number for wave members)
n_sele=25; %member equation approximation (the rest is stochastic)

tx=linspace(-x(end),x(end),2*length(x)-1); 

%functions preallocation with zeros
Fr=zeros(n_sele);
Am=zeros(n_sele);
y=zeros(length(x),n_sele);
fcorr=zeros(length(tx),n_sele);
ycorr=zeros(length(x),n_sele);

for n=1:n_sele
Fr(n)=x_peak(n); %Frequency
Am(n)=y_peak(n); %Amplitude
y_not_sync(:,n)=Am(n)*sin(2*pi*x*Fr(n)); %sin equations not in sync
fcorr(:,n)=xcorr(yv,y_not_sync(:,n)); %correlation function  
[y_fcorr(n),x_fcorr(n)]=max(fcorr(:,n)); %finding position of the correlation point
ycorr(:,n)=Am(n)*sin(2*pi*(x-tx(x_fcorr(n)))*Fr(n)); %sin equations in sync, tx(x_fcorr(n)) is the phase
ywave=sum(ycorr,2); %superimposition of single cycles
end

figure()
tiledlayout(2,1)
plot(nexttile,x,yv)
hold on
plot(nexttile,x,ywave)

SSres=sum((yv-ywave).^2,1);
SStot=sum((yv-mean(yv)).^2,1);
R2=1-SSres/SStot;

.

[William Rose]

Wow. @Star Strider found phase angles in the standard way, and it worked great, i.e. r^2=0.996.

@Giacomo, I now understand a little better than before what you did in your approach. You do the Fourier transform. Then you compute the amplitude spectrum and discard the phase information in the Fourier transform. You find 25 peaks in the amplitude spectrum. You use cross correlation to determine the integer sample offset that gives the best correlation between each sinusoidal component and the original signal. You use the cross-correlation results to apply a phase shift to each component used to recontruct the signal.

That is an interesting way to estimate phase angles. Why do you do it that way? I'm surte you have a reason. I never thought of determining phase angles that way. I think it is less accurate than using the phases determined directly from the Fourier transform, because your phase angles are limited to the discrete set of angles corresponding to integer numbers of samples.

You use a finite number of frequencies, smaller than all possible frequencies, to reconstruct a signal. Rather than use the first N harmonics, starting from the lowest frequency, you use harmonics which are peaks in the amplitude spectrum. That's a nice approach, because you recover a lot of signal energy with a limited number of frequency components.

The phase approach of @Star Strider is easier for me to understand, and it works very well.

[Giacomo]

@Star Strider thanks for your support. Now the interpolation works perfectly.

I need to ask you two more questions. One related to your code and one to a further implementation.

1- You calculated the phase as ph = angle(FFT). Why did you implement instead into the equation in this way?

for k = 1:numel(f)

Re(k,:) = A(k)*cos(2*pi*f(k)*x + angle(k));

end

2- Once i got the coefficients if i try to project into the future i got a flat profile after x(end)

x_proj=0:1:500;

for k=1:numel(f)

Re_proj(k,:)=A(k)*cos(2*pi*f(k)*x_proj+angle(k));

end

Res_proj=sum(Re_proj);

Is this approach not correct?

Thank you

[Star Strider]

As always, my pleasure!

1 The reconstructed parts of the Fourier transform of this signal are really strange. The sine and cosine terms appear to be completely out-of-phase, and while I would expect them to be separated by from each other, to have them completely out-of-phase () is not something I would expect. I generally do not do discrete since and cosine transforms, however according to the Wikipedia article (checked for reference) this analysis is correect. II changed my original approach to usee the sin and cos transforms instead of cos and phase, so the zero result is appropriate, although somewhat mystifyinig.

There is something here that I cannot explain, however thte new analysis works —

x =[0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19    20 21    22    23    24    25    26    27    28    29    30    31    32    33    34    35    36    37    38    39    40    41 42    43    44    45    46    47    48    49    50    51    52    53    54    55    56    57    58    59    60    61    62 63    64    65    66    67    68    69    70    71    72    73    74    75    76    77    78    79    80    81    82    83 84    85    86    87    88    89    90    91    92    93    94    95    96    97    98    99   100   101   102   103   104 105   106   107   108   109   110   111   112   113   114   115   116   117   118   119   120   121   122   123   124   125 126   127   128   129   130   131   132   133   134   135   136   137   138   139   140   141   142   143   144   145   146 147   148   149   150   151   152   153   154   155   156   157   158   159   160   161   162   163   164   165   166   167 168   169   170   171   172   173   174   175   176   177   178   179   180   181   182   183   184   185   186   187   188 189   190   191   192   193   194   195   196   197   198   199   200   201   202   203   204   205   206   207   208   209 210   211   212   213   214   215   216   217   218   219   220   221   222   223   224   225   226   227   228   229   230 231   232   233   234   235   236   237   238   239   240   241   242   243   244   245   246   247   248   249   250   251 252   253   254   255   256   257   258   259   260   261   262   263   264   265   266   267   268   269   270   271   272 273   274   275   276   277   278   279   280   281   282   283   284   285   286   287   288   289   290   291   292   293 294   295   296   297   298];
y =[-10.5109   -1.2195    7.8469   14.1429   16.2819   14.6567   11.2349    8.6033    8.7078   11.8968   16.7308   20.6538   21.1976   17.1386    9.0625   -0.9103  -10.0699  -16.2555  -18.6490  -17.8892  -15.5312  -13.1869  -11.8010  -11.3909  -11.2983  -10.7266   -9.2172   -6.8066   -3.8283   -0.5531    3.0506    7.2480   12.1846   17.4779   22.0492   24.3720   23.0792   17.6471    8.7990   -1.6092  -11.2030  -17.9301  -20.7384  -19.7912  -16.1598  -11.2007   -5.9639   -0.9310    3.8327    8.2740   11.9973   14.2207  14.1022   11.2560    6.1434    0.0715   -5.2581   -8.4492   -9.0052   -7.5011   -5.2318   -3.5181   -3.0353   -3.5248   -4.0378 -3.5722   -1.7391    0.9238    3.1897    3.8204    2.3237   -0.6352   -3.4286   -4.2130   -1.8732    3.3114    9.6375   14.6372    16.1367   13.2253    6.6788   -1.3565   -8.2684  -12.0471  -12.0361   -9.0487   -4.8460   -1.2646    0.5771    0.6378   -0.2435    -0.9455   -0.7268    0.3869    1.6692    2.1536    1.2313   -0.9304   -3.3716   -4.7455   -3.9615   -0.7664    4.0153    8.7097   11.4121   10.7500    6.4734   -0.3733   -7.7794  -13.4745  -15.7392  -14.0051   -9.0118   -2.4672    3.6402    7.7917    9.4438   9.0905    7.8818    6.9889    7.0208    7.7607    8.3378    7.7356    5.3731    1.4572   -3.0786   -7.0221   -9.4841  -10.3331 -10.2179  -10.1477  -10.8409  -12.1942  -13.1898  -12.3401   -8.4894   -1.5794    7.0370   14.9522   19.6282   19.4898   14.6770    7.0773   -0.4123   -5.0824   -5.5302   -2.1934    2.9356    7.2711    8.8090    6.9740    2.7630   -1.8432   -4.8836   -5.3432   -3.5700   -1.0217    0.4950   -0.3003   -3.5852   -8.3959  -13.0926  -16.0629  -16.3249  -13.7819   -9.0856   -3.2546    2.7213    8.1070   12.4336   15.3871   16.7195   16.2816   14.1588   10.7981    6.9895    3.6470    1.4554    0.5552    0.4452    0.1855  -1.1633   -4.0650   -8.2230  -12.6336  -15.9623  -17.0760  -15.4899  -11.5380   -6.2099   -0.7514    3.7744    6.7799    8.2032  8.3496    7.6510    6.4775    5.0813    3.6523    2.4035    1.5944    1.4562    2.0573    3.2000    4.4279    5.1684    4.9524 3.6069    1.3246   -1.4270   -4.1073   -6.3144   -7.8987   -8.9312   -9.5519   -9.7939   -9.4848   -8.2921   -5.8971   -2.2196  2.4193    7.2823   11.3738   13.7237   13.6958   11.2107    6.8034    1.4920   -3.5023   -7.0973   -8.6461   -8.0951   -5.9509 -3.0666   -0.3190    1.7049    2.8917    3.5746    4.2977    5.4674    7.0454    8.4417    8.6881    6.8505    2.5180   -3.8545 -10.9311  -16.8617  -19.9141  -19.1236  -14.6870   -7.9182   -0.7697    4.8922    7.9559    8.4032    7.1844    5.7055    5.1666   6.0469    7.9541    9.8788   10.7135    9.7826    7.1424    3.5299    0.0083   -2.5093   -3.6337   -3.6112   -3.1454   -3.0333 -3.7945   -5.4569   -7.5774   -9.4631  -10.4763  -10.2829   -8.9475   -6.8554   -4.5183   -2.3599   -0.5720    0.9113    2.3268 3.9330    5.8640    8.0423   10.1782   11.8523   12.6467   12.2813   10.7082    8.1341    4.9623    1.6760   -1.2969   -3.6892];
% yv=y';
% Length=x(end);
% Datapoints=length(x);
% fs=round(Datapoints/Length);

fs = 1/(x(2) - x(1))

fs = 1

L = numel(x)

L = 299

nfft = 2^nextpow2(L)*2;

FFT = fft(y(:),nfft)/L;
% FFT = 2*abs(FFT(1:nfft/2+1));
f = fs/2*linspace(0,1,nfft/2+1);

A = real(FFT(1:numel(f)));
B = imag(FFT(1:numel(f)));
ph = angle(FFT(1:numel(f)));

figure
plot(f, A, 'DisplayName','A (Amplitude units)')
hold on
plot(f, B, 'DisplayName','B (Amplitude units)')
plot(f, ph, 'DisplayName','\Phi (radians)')
hold off
grid
xlim([0 0.15])
legend('Location','best')
xlabel('Frequency')
ylabel('Component')
title('Coefficients & Phase vs. Frequency')

Re = A(:).*cos(2*pi*f(:)*x);
Im = B(:).*sin(2*pi*f(:)*x);

figure
plot(x, sum(Re), 'DisplayName','Re')
hold on
plot(x,sum(Im), 'DisplayName','Im')
plot(x, sum(Re)+sum(Im), 'DisplayName','Re+Im')
hold off
grid
legend('Location','best')

Res = (sum(Re)-sum(Im))/2

Res = 1×299

-8.9085 -0.9369 6.8093 12.2162 14.0312 12.6561 9.7100 7.4732 7.5463 10.2931 14.4155 17.7907 18.2400 14.7811 7.8500 -0.6722 -8.5309 -13.8106 -15.8762 -15.2093 -13.2068 -11.1833 -10.0130 -9.6455 -9.5826 -9.0768 -7.8008 -5.7205 -3.1869 -0.3664

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Mdl = fitlm(y, Res)

Mdl = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate        SE        tStat    pValue
                   ________    __________    _____    ______

    (Intercept)    0.098988    0.00047522    208.3      0   
    x1              0.85619    5.1073e-05    16764      0   


Number of observations: 299, Error degrees of freedom: 297
Root Mean Squared Error: 0.00821
R-squared: 1,  Adjusted R-Squared: 1
F-statistic vs. constant model: 2.81e+08, p-value = 0

figure
plot(x, Res, 'LineWidth',1.5, 'DisplayName','Reconstructed Signal')
hold on
plot(x, y, '--', 'LineWidth',1, 'DisplayName','Original Signal')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

[K,X] = ndgrid(1:numel(f),x);

figure
plot3(X, K, Re.')
grid on
xlabel('Time')
ylabel('Term')
title('Re')

return

[y_peak,x_peak]=findpeaks(FFT,f);

%wave equation yn=sum(An*sin(2*pi*Fn*(x+Pn)) n=number of peaks (number for wave members)
n_sele=25; %member equation approximation (the rest is stochastic)

tx=linspace(-x(end),x(end),2*length(x)-1); 

%functions preallocation with zeros
Fr=zeros(n_sele);
Am=zeros(n_sele);
y=zeros(length(x),n_sele);
fcorr=zeros(length(tx),n_sele);
ycorr=zeros(length(x),n_sele);

for n=1:n_sele
Fr(n)=x_peak(n); %Frequency
Am(n)=y_peak(n); %Amplitude
y_not_sync(:,n)=Am(n)*sin(2*pi*x*Fr(n)); %sin equations not in sync
fcorr(:,n)=xcorr(yv,y_not_sync(:,n)); %correlation function  
[y_fcorr(n),x_fcorr(n)]=max(fcorr(:,n)); %finding position of the correlation point
ycorr(:,n)=Am(n)*sin(2*pi*(x-tx(x_fcorr(n)))*Fr(n)); %sin equations in sync, tx(x_fcorr(n)) is the phase
ywave=sum(ycorr,2); %superimposition of single cycles
end

figure()
tiledlayout(2,1)
plot(nexttile,x,yv)
hold on
plot(nexttile,x,ywave)

SSres=sum((yv-ywave).^2,1);
SStot=sum((yv-mean(yv)).^2,1);
R2=1-SSres/SStot;

That is the best I can do with this.

.

[Giacomo]

I'm sorry but the problem remain the same.

Even with these calculations the function still shows a flat profile above x(end)=298

x1=0:1:500;

Re1 = A(:).*cos(2*pi*f(:)*x1);

Im1 = B(:).*sin(2*pi*f(:)*x1);

Res1 = (sum(Re1)-sum(Im1))/2;

figure(1)

plot(x1,Res1)

The interpolation looks great but the projection doesn't work.

Thanks

[Star Strider]

This is interesting, although not enlightening.

When I plotted all the data, all the elements of the Fourier series are eevaluated, apparently correctly. Summing them leads to interesting behaviour that I was previously not aware of. After the end of the original ‘x’ vector, the summed elements apparently cancel, giiving a decaying sinusoiidal oscillation to zero. The more elements there are, the more quickly the sum of the vectors reaches zero.

There has to be an analytic explanation for this, however it currently eludes me. It would seem that the extended time vector (‘x2’ here) would endlessly repeat the original series for all time, however that obviously does not occur. All the elements of the Fourier series continue to be present, however they sum to zero beyond the original independent variable range.

My analysis —

x =[0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19    20 21    22    23    24    25    26    27    28    29    30    31    32    33    34    35    36    37    38    39    40    41 42    43    44    45    46    47    48    49    50    51    52    53    54    55    56    57    58    59    60    61    62 63    64    65    66    67    68    69    70    71    72    73    74    75    76    77    78    79    80    81    82    83 84    85    86    87    88    89    90    91    92    93    94    95    96    97    98    99   100   101   102   103   104 105   106   107   108   109   110   111   112   113   114   115   116   117   118   119   120   121   122   123   124   125 126   127   128   129   130   131   132   133   134   135   136   137   138   139   140   141   142   143   144   145   146 147   148   149   150   151   152   153   154   155   156   157   158   159   160   161   162   163   164   165   166   167 168   169   170   171   172   173   174   175   176   177   178   179   180   181   182   183   184   185   186   187   188 189   190   191   192   193   194   195   196   197   198   199   200   201   202   203   204   205   206   207   208   209 210   211   212   213   214   215   216   217   218   219   220   221   222   223   224   225   226   227   228   229   230 231   232   233   234   235   236   237   238   239   240   241   242   243   244   245   246   247   248   249   250   251 252   253   254   255   256   257   258   259   260   261   262   263   264   265   266   267   268   269   270   271   272 273   274   275   276   277   278   279   280   281   282   283   284   285   286   287   288   289   290   291   292   293 294   295   296   297   298];
y =[-10.5109   -1.2195    7.8469   14.1429   16.2819   14.6567   11.2349    8.6033    8.7078   11.8968   16.7308   20.6538   21.1976   17.1386    9.0625   -0.9103  -10.0699  -16.2555  -18.6490  -17.8892  -15.5312  -13.1869  -11.8010  -11.3909  -11.2983  -10.7266   -9.2172   -6.8066   -3.8283   -0.5531    3.0506    7.2480   12.1846   17.4779   22.0492   24.3720   23.0792   17.6471    8.7990   -1.6092  -11.2030  -17.9301  -20.7384  -19.7912  -16.1598  -11.2007   -5.9639   -0.9310    3.8327    8.2740   11.9973   14.2207  14.1022   11.2560    6.1434    0.0715   -5.2581   -8.4492   -9.0052   -7.5011   -5.2318   -3.5181   -3.0353   -3.5248   -4.0378 -3.5722   -1.7391    0.9238    3.1897    3.8204    2.3237   -0.6352   -3.4286   -4.2130   -1.8732    3.3114    9.6375   14.6372    16.1367   13.2253    6.6788   -1.3565   -8.2684  -12.0471  -12.0361   -9.0487   -4.8460   -1.2646    0.5771    0.6378   -0.2435    -0.9455   -0.7268    0.3869    1.6692    2.1536    1.2313   -0.9304   -3.3716   -4.7455   -3.9615   -0.7664    4.0153    8.7097   11.4121   10.7500    6.4734   -0.3733   -7.7794  -13.4745  -15.7392  -14.0051   -9.0118   -2.4672    3.6402    7.7917    9.4438   9.0905    7.8818    6.9889    7.0208    7.7607    8.3378    7.7356    5.3731    1.4572   -3.0786   -7.0221   -9.4841  -10.3331 -10.2179  -10.1477  -10.8409  -12.1942  -13.1898  -12.3401   -8.4894   -1.5794    7.0370   14.9522   19.6282   19.4898   14.6770    7.0773   -0.4123   -5.0824   -5.5302   -2.1934    2.9356    7.2711    8.8090    6.9740    2.7630   -1.8432   -4.8836   -5.3432   -3.5700   -1.0217    0.4950   -0.3003   -3.5852   -8.3959  -13.0926  -16.0629  -16.3249  -13.7819   -9.0856   -3.2546    2.7213    8.1070   12.4336   15.3871   16.7195   16.2816   14.1588   10.7981    6.9895    3.6470    1.4554    0.5552    0.4452    0.1855  -1.1633   -4.0650   -8.2230  -12.6336  -15.9623  -17.0760  -15.4899  -11.5380   -6.2099   -0.7514    3.7744    6.7799    8.2032  8.3496    7.6510    6.4775    5.0813    3.6523    2.4035    1.5944    1.4562    2.0573    3.2000    4.4279    5.1684    4.9524 3.6069    1.3246   -1.4270   -4.1073   -6.3144   -7.8987   -8.9312   -9.5519   -9.7939   -9.4848   -8.2921   -5.8971   -2.2196  2.4193    7.2823   11.3738   13.7237   13.6958   11.2107    6.8034    1.4920   -3.5023   -7.0973   -8.6461   -8.0951   -5.9509 -3.0666   -0.3190    1.7049    2.8917    3.5746    4.2977    5.4674    7.0454    8.4417    8.6881    6.8505    2.5180   -3.8545 -10.9311  -16.8617  -19.9141  -19.1236  -14.6870   -7.9182   -0.7697    4.8922    7.9559    8.4032    7.1844    5.7055    5.1666   6.0469    7.9541    9.8788   10.7135    9.7826    7.1424    3.5299    0.0083   -2.5093   -3.6337   -3.6112   -3.1454   -3.0333 -3.7945   -5.4569   -7.5774   -9.4631  -10.4763  -10.2829   -8.9475   -6.8554   -4.5183   -2.3599   -0.5720    0.9113    2.3268 3.9330    5.8640    8.0423   10.1782   11.8523   12.6467   12.2813   10.7082    8.1341    4.9623    1.6760   -1.2969   -3.6892];
% yv=y';
% Length=x(end);
% Datapoints=length(x);
% fs=round(Datapoints/Length);

fs = 1/(x(2) - x(1))

fs = 1

L = numel(x)

L = 299

nfft = 2^nextpow2(L)*2;

FFT = fft(y(:),nfft)/L;
% FFT = 2*abs(FFT(1:nfft/2+1));
f = fs/2*linspace(0,1,nfft/2+1);

A = real(FFT(1:numel(f)));
B = imag(FFT(1:numel(f)));
ph = angle(FFT(1:numel(f)));

figure
plot(f, A, 'DisplayName','A (Amplitude units)')
hold on
plot(f, B, 'DisplayName','B (Amplitude units)')
plot(f, ph, 'DisplayName','\Phi (radians)')
hold off
grid
xlim([0 0.15])
legend('Location','best')
xlabel('Frequency')
ylabel('Component')
title('Coefficients & Phase vs. Frequency')

Re = A(:).*cos(2*pi*f(:)*x);
Im = B(:).*sin(2*pi*f(:)*x);


figure
plot(x, sum(Re), 'DisplayName','Re')
hold on
plot(x,sum(Im), 'DisplayName','Im')
plot(x, sum(Re)+sum(Im), 'DisplayName','Re+Im')
hold off
grid
legend('Location','best')

Res = (sum(Re)-sum(Im))/2

Res = 1×299

-8.9085 -0.9369 6.8093 12.2162 14.0312 12.6561 9.7100 7.4732 7.5463 10.2931 14.4155 17.7907 18.2400 14.7811 7.8500 -0.6722 -8.5309 -13.8106 -15.8762 -15.2093 -13.2068 -11.1833 -10.0130 -9.6455 -9.5826 -9.0768 -7.8008 -5.7205 -3.1869 -0.3664

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Mdl = fitlm(y, Res)

Mdl = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate        SE        tStat    pValue
                   ________    __________    _____    ______

    (Intercept)    0.098988    0.00047522    208.3      0   
    x1              0.85619    5.1073e-05    16764      0   


Number of observations: 299, Error degrees of freedom: 297
Root Mean Squared Error: 0.00821
R-squared: 1,  Adjusted R-Squared: 1
F-statistic vs. constant model: 2.81e+08, p-value = 0

figure
plot(x, Res, 'LineWidth',1.5, 'DisplayName','Reconstructed Signal')
hold on
plot(x, y, '--', 'LineWidth',1, 'DisplayName','Original Signal')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

% [K,X] = ndgrid(1:numel(f),x);
% 
% figure
% plot3(X, K, Re.')
% grid on
% xlabel('Time')
% ylabel('Term')
% title('Re')

x2 = 0:500;
Re2 = A(:).*cos(2*pi*f(:)*x2);
Im2 = B(:).*sin(2*pi*f(:)*x2);

Q1 = A(2)*cos(2*pi*f(2)*x2);

% Q2 = Re2(1:10,:)

% Qir = reshape(Q1(1:500), 10, [])

figure
% plot(x2, Re2(1:10,:))
plot(x2, Re2)
grid

% ixr = 295:375;
ixr = 291:350;

figure
% plot(x2, Re2(1:10,:))
plot(x2(ixr), Re2(2:106,ixr))
grid

figure
% plot(x2, Re2(1:10,:))
plot(x2(ixr), sum(Re2(2:106,ixr))  )
grid

figure
% plot(x2, Re2(1:10,:))
plot(x2, sum(Re2))
grid

return

Res2 = (sum(Re2)-sum(Im2))/2
Res2 = sum(Re2)

figure
plot(x2, Res2, 'LineWidth',1.5, 'DisplayName','Reconstructed Signal')
hold on
% plot(x, y, '--', 'LineWidth',1, 'DisplayName','Original Signal')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

return

[y_peak,x_peak]=findpeaks(FFT,f);

%wave equation yn=sum(An*sin(2*pi*Fn*(x+Pn)) n=number of peaks (number for wave members)
n_sele=25; %member equation approximation (the rest is stochastic)

tx=linspace(-x(end),x(end),2*length(x)-1); 

%functions preallocation with zeros
Fr=zeros(n_sele);
Am=zeros(n_sele);
y=zeros(length(x),n_sele);
fcorr=zeros(length(tx),n_sele);
ycorr=zeros(length(x),n_sele);

for n=1:n_sele
Fr(n)=x_peak(n); %Frequency
Am(n)=y_peak(n); %Amplitude
y_not_sync(:,n)=Am(n)*sin(2*pi*x*Fr(n)); %sin equations not in sync
fcorr(:,n)=xcorr(yv,y_not_sync(:,n)); %correlation function  
[y_fcorr(n),x_fcorr(n)]=max(fcorr(:,n)); %finding position of the correlation point
ycorr(:,n)=Am(n)*sin(2*pi*(x-tx(x_fcorr(n)))*Fr(n)); %sin equations in sync, tx(x_fcorr(n)) is the phase
ywave=sum(ycorr,2); %superimposition of single cycles
end

figure()
tiledlayout(2,1)
plot(nexttile,x,yv)
hold on
plot(nexttile,x,ywave)

SSres=sum((yv-ywave).^2,1);
SStot=sum((yv-mean(yv)).^2,1);
R2=1-SSres/SStot;

.

[Giacomo]

Hi @William Rose Me and @Star Strider are facing a major issue. We succesfully interpolated the data according to this Fourier Model f(x) =sum[ A(k)*cos(2*pi*f(k)*x + angle(k))]. However when we try to project the function further in future (having already calculated model coefficients) the function tends to an high frequency-small amplitude function (stochastic beahaviour, like we plot only the background noise without the wave meat). Do you have any suggestion?

I attach the model.

Thanks

[William Rose]

@Giacomo,

The results you are observing are due to a combination of factors:

• you use nfft = 2^nextpow2(L)*2;
• you reconstruct the signal using only the real parts of the fft

By using nfft=1024, when the signal has length 299, the signal is zero-padded out to 1024 samples. Thus the inverse FFT, if you computed it, would give a periodic signal with period 1024 and about 725 zeros. But you do not do a standard inverse FFT. You reconstruct the signal using only the real parts, which alters the reconstructed signal in other ways: 1. the reconstruction (whose period is 1024) is imperfect, and 2. the reconstruction is symmetric about x=zero, and (equivalently) is reflected about 511.5 (i.e. points 512-1023 are a backwards version of points 0 to 511).

I have attached 2 scripts. The first script is fourierModel2.m. This is the model you attached most recently, with 2 changes: I added clear at the start, and I changed x1 from 0:500 to 0:1500. By increasing x1 this way, you can see the symmetry and periodic nature of the reconstrcucted signal. The figure below is from fourierModel2.m.

The second script is fourierModel3.m. Here I use nfft=L=299, and I do a standard reconstruction of the signal, using the complex coefficients of the FFT, instead of using the real parts only. The plot below shows the original signal, the reconstructed signal using real part only, and the reconstructed signal using the complex Fourier coefficients. Note that the complex reconstruction is exactly correct, the real reconstru ction is not very accurate, and the real and the complex reconstruction are periodic with period 299.

[William Rose]

@Giacomo, I noticed a mistake in the script Fourier model.txt, which you attached to your recent posting. The angle(k) term is wrong.

A = real(FFT(1:numel(f)));
ph = angle(FFT);
for k = 1:numel(f)
    Re(k,:) = A(k)*cos(2*pi*f(k)*x + angle(k));
end

The corrected code (using real parts, as you do) would be

A = real(FFT(1:numel(f)));
ph = angle(FFT);
for k = 1:numel(f)
    Re(k,:) = A(k)*cos(2*pi*f(k)*x + ph(k));
end

The script fourierModel2.m, which I posted above, was mostly a copy of your script, so it included the same error. When I correct the error in fourierModel2.m (see attached script fourierModel2a.m), the signal reconstructed with the real parts fits the data much less well than before. I don't know why the reconstruction with real parts only (which is wrong) combined with the wrong phase angles resulted in a better fit than a reconstruction with real parts plus the correct phase angles.

If you use the magnitudes plus the correct phase angles, you get an exact reconstruction of the signal. See attached script 3a which uses

Ar = real(X);
Am = abs(X);
ph = angle(X);
x1=0:500;
for k = 1:nfft
    Re(k,:) = Ar(k)*cos(2*pi*f(k)*x1 + ph(k));
    Mg(k,:) = Am(k)*cos(2*pi*f(k)*x1 + ph(k));
end
Res = sum(Re);
Mag = sum(Mg);

[I'm having trouble arttaching code and inserting screenshot. I hope the explanation above is sufficient even with out the full script and screenshot.]

Answer 2

[William Rose]

@Giacomo,

I do not understand everything you are doing in your script. But it looks lto me like you discard the phase information. I think that if you retain and use the phase information, you will get better results.

By the way you wrote, in th comments of your script,

%wave equation yn=sum(An*sin(2*pi*Fn*(x+Pn)) n=number of peaks (number for wave members)

I would call that the equation for Fourier synthesis. The wave equation, to me and to many others, is

or a 2D or 3D version of the above equation.

[Giacomo]

I implemented the phase with the correlation function. Basically the script calculates An*sin(2*pi*Fn*x) then i use the autocorrelation with the original data to find the phase. My new sync equation will be An*sin(2*pi*Fn*(x-delay)).

A wave equation is the sum of all its cyclicals. I dont need a 2d or 3d model. The superimposition of all cycles is just fine.