Transform disharmonic signal to harmonic

Hello Matlab community,

I'd like to convert an aperiodic disharmonic measurement signal into a periodic harmonic signal. The signal given is a distance(x) over displacement(z) signal. Via hanning I've converted it into a periodic signal. Now for the disharmonic part: The idea is to split the signal into n small pieces and to reproduce it via fft(z) into a sine signal, containing data of frequency, amplitude and distance of the event, in order to get a complete sine function for every single of the n pieces I've split the original signal into. This is all for the purpose of solving the differential equation of a mass oscillator with harmonic force (or displacement) input. In the end I'd like to get a converted measurement signal consisting of n different sine functions, where I could store the data of frequency, amplitude and phase in an array, well at least that's the idea... I'd be very thankful for any idea as to how to generate the fft out of this. Will I need a 3-dimensional array or could this be solved with angle(z) for the distance?

Answer 1

[Star Strider]

I am a bit confused, so I may need you to clarify your question.

It seems to me that if you want to create a ‘long’ version of a ‘short’ sequence, you need to get the amplitudes and frequencies of the real (cosine) terms and complex (sine) terms of your ‘short’ sequence with the Fourier transform. Then create a time vector for as long a time as you want, and give the time vector as the argument to your sum-of-sines-and-cosines function to create the signal you want.

[Tobias Engelhardt]

Hello Star Strider, thank you for your help.

I have given a long signal. So I look at it in small parts where each part may be described as a sine function. To fit this sine function as precise as possible to this small part of my signal, the amplitude and the frequency is needed. After that, I'd like to fit all the small functions I've generated back together, to rebuild my input signal in sine functions. Therefor I think I'd need a timestamp or some mark as to where the viewed element n belongs to. Come to think that I could do this in a loop:

for(k=1:1:length(distance_x))
   freq(k)=fft(displacement_z(k));   %get frequency data of signal z at point k
   ampl(k)=abs(freq);   %get amplitude data of signal z at point k
end
input=ampl*(cos(freq*time));  %I think cos() is better for fitting later

To get the time vector I'd need to set a speed v, but that is a minor problem. Not sure about the amplitude data though.

[Star Strider]

My pleasure.

I’m still not sure I understand what you want to do.

You cannot get any meaningful information from the Fourier transform of a single point, since the Fourier transform of a single point is the single point itself.

See the documentation for the fft (link) function for details. (I like the R2015a version documentation, so I linked to it here.)

The Fourier transform will isolate the amplitudes at each frequency, so if a particular frequency is element freq(k) from the frequency vector created in your code (as described in the documentation), and the Fourier transform coefficients at that frequency are A(k)+j*B(k) (where ‘j’ is the imaginary operator and ‘A’ and ‘B’ are real numbers) and ‘t’ is your vector of times, the time-domain curve corresponding to that frequency is:

y = @(t,k) A(k)*cos(freq(k)*t) + B(k)*sin(freq(k)*t);               % For Frequency In ‘radians/time’
y = @(t,k) A(k)*cos(2*pi*freq(k)*t) + B(k)*sin(2*pi*freq(k)*t);     % For Frequency In ‘cycles/time’

That is how I would synthesise them.

[Tobias Engelhardt]

I'm sorry to ask again for your help, but I could not get to a solution yet. There must be a function to get an interpolated harmonic signal for my inharmonic input signal, like this:

The blue disharmonic signal y=[0 3 7 6 1 0] can be displayed closely as the function f=7*sin(0.1*2*pi*x) Therefor I want to extract the frequency of the inharmonic signal over let's say 5 steps and get the amplitude along with it. Which would be freq=0.1 and apml=7

Maybe you've already posted the correct code and I just don't get it... In that case I'm sincerely sorry.

[Star Strider]

If you want a least-squares ‘best fit’ of a sine function to your data, see if this does what you want:

x = 0:5;
y = [0 3 7 6 1 0];
f = @(b,x) b(1).*sin(b(2)*2*pi*x);
normrsd = @(b) norm(y - f(b,x));
[B, nrmrsd] = fminsearch(normrsd, [7; 0.1]);
fprintf(1, '\n\tAmplitude =\t%5.2f\n\tFrequency =\t%5.2f\n\n', B)

xplot = linspace(min(x), max(x));

figure(1)
plot(x, y)
hold on
plot(xplot, f(B,xplot))
hold off
grid

  Amplitude =   6.07
  Frequency =   0.11

[Tobias Engelhardt]

You Sir, are a genius!

[Star Strider]

Thank you!

[Tobias Engelhardt]

Just one more little thing: Since I want to do this operation on a stochastic signal which gets split up in little parts, it would be great if the interpolated sine function would end at the same amplitude as the little piece of my signal, so I can put the signal back together with a smooth connection between the parts. How could this be achieved?

Thank you.

[Star Strider]

It amazes me what it’s possible to do with the fmincon function (probably because I don’t use it often, so don’t have much experience with it).

I was quite happy with the result it gave here:

x = 0:5;
y = [0 3 7 6 1 0];
f = @(b,x) b(1).*sin(b(2)*2*pi*x);                                      % Objective Function
normrsd = @(b) norm(y - f(b,x));                                        % Norm-Of-Residuals Function (To Be Optimised)
% % [B, nrmrsd] = fminsearch(normrsd, [7; 0.1]);
c = @(b) [y(end) - sin(b(2)*2*pi*x), -b(1)];                            % Inequality Constraints — The Function Has To Be ≥ ‘y(end)’
ceq = @(b)  [];                                                         % Equality Constraints (None)
b2con = @(b) deal(c(b),ceq(b));                                         % Create Constraint Function
[B,nrmrsd] = fmincon(normrsd, [7; 0.1], [],[],[],[],[],[],b2con);       % Estimate Parameters Given Constraints
fprintf(1, '\n\tAmplitude =\t%5.2f\n\tFrequency =\t%5.2f\n\n', B)

xplot = linspace(min(x), max(x));
fit_f = f(B,xplot);

figure(1)
plot(x, y)
hold on
plot(xplot, f(B,xplot))
hold off
grid
xlabel('\itx\rm \rightarrow')
ylabel('\ity\rm \rightarrow')
legend('Original Signal', 'Sine Function Approximation', 'Location','S')

I documented the code this time with comments.

The nonlinear constraint function construction needs a more detailed description, since it’s not obvious.

The inequality constraint function ‘c(b)’ derives from:

sin(b(2)*2*pi*x) ≥ y(end)

that becomes:

y(end) - sin(b(2)*2*pi*x)

and since (here) b(1) simply has to be ≥ 0, it becomes simply:

-b(1)

See the documentation for Nonlinear Constraints for details on the conventions for creating the constraint functions.

My code produces ...

  Amplitude =   5.89
  Frequency =   0.10

The Plot —