Plot Fourier Transform (FFT) Interpolation Model into future

2 ビュー (過去 30 日間)
Giacomo
Giacomo 2024 年 10 月 1 日
コメント済み: Giacomo 2024 年 10 月 3 日
Hi i'm gonna open another discussion hoping more people will join it and someone will find the solution. I have a set of data which have been interpolated through Fourier Transform to build a function f(x)=A*cos(2*pi*F*x+P). All the coefficients have been succesfully calculated and the function fit perfectly the data. The problem is that when i try to plot the equation to further x value the function tends immediately to a stochastic behaviour (small amplitude-high frequency wave) which it doesn't make any sense.
In other words after the interpolation period the function collapse immediately into its noise and strong amplitudes disappear.
Does anyone can solve this problem?
I attach the model.
Thanks for you all in advance
  1 件のコメント
Star Strider
Star Strider 2024 年 10 月 1 日
The Fourier coeefficients create extrapolated values appropriate to an independent variable vector that is a direct extension (0 to 500) of the original (0 to 299). Beyond the range of original independent variable vector, the created waves all appear to completely cancel, creatiing a constant line instead of extrapolatiing the original siignal.
A theoretical discussion of that behaviour would be appreciiated, since the reason is not obvious (at least to me).

サインインしてコメントする。

回答 (2 件)

Mathieu NOE
Mathieu NOE 2024 年 10 月 1 日
編集済み: Mathieu NOE 2024 年 10 月 1 日
hello Giacomo
first I corrected some minor bugs (before we dig into more technical stuff)
  • A is an amplitude so you should take abs(fft) and not real(fft)
  • the phase you use is wrong in Re(k,:) = A(k)*cos(2*pi*f(k)*x + angle(k)) must be replaced by Re1(k,:) = A(k)*cos(2*pi*f(k)*x1 + ph(k)) , as you have computed a few lines above : ph = angle(FFT);
  • fft amplitude normalisation; there is always lot of discussion of how one should normalize the fft output to get the right amplitude; for me you should replace FFT = fft(y(:),nfft)/L; with : FFT = fft(y(:),nfft)/(nfft/2); I have to dig again into this topic to give the correct info , but I am a bit in a hurry today...
  • next point I do not get is why you have intentionnally increased the nfft size (the factor 2 at the trailing end ) : nfft = 2^nextpow2(L)*2; the code works fine without it , so below we have : nfft = 2^nextpow2(L) = 512 here;
Now , we have something that works , or seem to work as we have correctly reconstructed the original signal. Unfortunately the signal original length is only 299 samples long , and with the above corrections (nfft = 2^nextpow2(L), nfft = 512
the thing is , when you do the signal generation, there will be almost zero output (almost because there is no "perfect" reconstruction here ) for samples between 299+1 and 512-1 , then the signal is again present (a simple 512 sampls shifted version of the first one)
if you cerate enough samples you would see this phenomenon
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];
fs = 1/(x(2) - x(1));
L = numel(x);
nfft = 2^nextpow2(L);
FFT = fft(y(:),nfft)/(nfft/2);
f = fs/2*linspace(0,1,nfft/2+1);
A = abs(FFT(1:numel(f)));
ph = angle(FFT(1:numel(f)));
for k = 1:numel(f)
Re(k,:) = A(k)*cos(2*pi*f(k)*x + ph(k));
end
Res = sum(Re);
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')
x1=0:1:2000;
for k = 1:numel(f)
Re1(k,:) = A(k)*cos(2*pi*f(k)*x1 + ph(k));
end
Res1=sum(Re1);
figure()
plot(x1,Res1)
  10 件のコメント
Mathieu NOE
Mathieu NOE 2024 年 10 月 3 日
hello
maybe I didn't catch the goal here - all I can say is that depending on your model choice you get what you paid for... if you use Fourier series , you can only expect to get correct results and predictions for cyclic and stationnary signals.
If you need a more sophisticated approach to predict non periodic / unstationnary signals, I would rather consider ARMAX model or ANN / Deap Learning
but as I onced heard (or read) , prediction is a difficult job , especially about the future :)
so even the best model cannot predict the future with high accuracy forever.
Giacomo
Giacomo 2024 年 10 月 3 日
Thank you Paul for this last consideration about periodicity. Didn't notice it. It is not a problem.
Both of you have been of great help.
Thanks
Giacomo

サインインしてコメントする。


Paul
Paul 2024 年 10 月 1 日
編集済み: Paul 2024 年 10 月 1 日
I've only very briefly skimmed this code and the other thread. My guess is that you're seeing the effect of zero-padding the FFT when you go back to the time domain.
Let's take a simple example of a 10-point sequence and zero pad to 20 points on the call to fft
rng(100);
x = rand(1,10);
X = fft(x,20);
Going back to the time domain and plotting/comparing
xr = ifft(X,'symmetric');
figure
hold on
stem(0:9,x,'x')
stem(0:19,xr)
legend('x','xr')
The result for xr is exactly as expected based on the zero-padding the input on the call to fft(), which is equivalent to
% fft([x zeros(1,10)]);
If you were to implement the equation for xr explicitly as a function of X (instead of using ifft() ) and use it to compute xr for n > 19, you sould see the periodic extension of xr, and the period will be N = 20.
So, in your case, see what happens if you don't zero pad on the call to fft.
Also, if you remove the zero padding to the nextpow2 you'll have to verify that the rest of the code is correct if the original signal has an odd number of elements (assuming odd-length inputs are of interest).
If, as seems to be the case, the only objective is to generate the periodic extension of the input, there are easier ways to do this than using fft() and developing your own series summation. OTOH, what you're doing is a good way to develop understanding of the underlying concepts.
p.s. IMO the code in "Fourier model.txt" is short enough that it would be better to just copy/paste it into your question (please use the code formatting (left most incone in the Code section in the toostrip) so it's easier for others to see and work with.
  1 件のコメント
Giacomo
Giacomo 2024 年 10 月 1 日
Hi Paul, thanks for your answer and to share with me all knowledge. I listened to your suggestion and answered to Mathieu (above) integrating Matlab code. Hope that it is more clear for you understanding my real problem and my final purpose. I would really appreciate if you could give it a look.
Thanks
Giacomo

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeVibration Analysis についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by