Improving personal Fourier Interpolation Model

6 次查看（过去 30 天）
Giacomo
Giacomo 2024-9-29
评论： William Rose 2024-10-2
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

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2024-9-29
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.
.
  10 个评论
显示 8更早的评论
William Rose
William Rose 2024-10-1
@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
William Rose 2024-10-2
@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.]

请先登录，再进行评论。

更多回答（1 个）

William Rose
William Rose 2024-9-29
@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.
  1 个评论
Giacomo
Giacomo 2024-9-29
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by