Construct phase spectrum before IFFT

Dear all,

I would like to produce a signal in time domain from a given set of regular waves in the frequency domain. Say 4 regular waves with known wave frequency, amplitude and phase.

To understand how the commands work, I plotted all four waves in a time sequence, used FFT on this time sequence and plotted the amplitude (abs(FFT_output)*2/length(time_series)) and phase spectrum (angle(FFT_output)*2/length(time_series)), both scaled by 2/length(time_series). Building up the amplitude spectrum from my known frequencies and amplitudes in the frequency domain, I insert the amplitudes in a vector with the same length as the frequency vector, at the positions of the corresponding frequency and use fliplr to make the spectrum symmetric. This works fine and corresponds exactly to the outputted amplitude spectrum from the FFT command.

For the phase spectrum I perform the same trick, insert all known phases into a vector at the positions of the frequency they belong to. This constructed phase spectrum however does in no way correspond to the phase spectrum constructed from the FFT output. The phase spectrum outputted by the FFT command contains phase information for all frequencies in the spectrum, though it is only an addition of 4 waves. How should I interpret this?

[David Goodmanson]

Hi Roderick, could you provide an example of your code including the time domain input?

One thing going on for sure is that you don't want to divide the angle by 2*length(...) but your second-to-last sentence means that the nature of the input signal is of interest.

[Roderick]

Hi David,

Thanks for your response. Please see the code below.

%% Signal definition omega=[1 2 4 8]; amplitude=[4 5 6 18]; phase=([1.7 0.3 1.15 0.6])*(pi);

Max_t=10; dt=0.01; t=(0:dt:Max_t); % Time vector df=1/(t(end)-t(1)); f=(0:length(t)-1)*df; % Frequency vector

%% Adding regular waves for time signal manually time_signal=zeros(length(omega),length(t)); for i=1:length(omega) time_signal(i,:)=amplitude(i)*cos((omega(i)*2*pi)*t+(phase(i))); end

T_addition=sum(time_signal,1);

%% Fast Fourier Transform of time signal Y=fft(T_addition);

ai=abs((Y*2)/(length(T_addition))); % Amplitude spectrum generated by FFT ph=angle((Y*2)/length(T_addition)); % Phase spectrum generated by FFT

% Re-construct Y when total amplitude spectrum is known Ya=(((ai.*(exp(1i*ph))).*(length(ai))/2));

T_fft_ifft=ifft(Ya);

%% Constructing amplitude and phase spectrum from known parameters error=0.01;

index=(1:length(omega))';

for i=1:length(omega) index(i,1)=find(f > (omega(i)-error) & f < (omega(i)+error)); end

Y_ai_parameters=zeros(length(t)); Y_ph_parameters=zeros(length(t)); scaling_factor=length(t)/2;

for i=1:length(index) Y_ai_parameters(index(i),1)=amplitude(i); Y_ai_parameters((end-index(i)-2),1)=amplitude(i); Y_ph_parameters(index(i),1)=phase(i); Y_ph_parameters((end-index(i)-2),1)=phase(i); end

Y_parameters=(((Y_ai_parameters.*(exp(1i*Y_ph_parameters)))));

T_parameters=ifft(Y_parameters);

%% Plotting figure subplot(4,1,1) plot(f,ai); title('Amplitude spectrum from FFT'); ylabel('ai'); subplot(4,1,2) plot(f,ph); title('Phase spectrum from FFT'); ylabel('ph'); subplot(4,1,3) plot(f,Y_ai_parameters, '-r') hold on %plot(f_half,ai_half) title('Amplitude spectrum from parameters'); subplot(4,1,4) plot(f,Y_ph_parameters, '-r') hold on %plot(f_half,ph_half) title('Phase spectrum from parameters');

figure subplot(3,1,1) plot(t,T_addition); title('Time signal by adding regular waves'); xlabel('T [s]'); ylabel('A'); subplot(3,1,2) plot(t,T_fft_ifft); title('Time signal from IFFT'); xlabel('T [s]'); ylabel('A'); subplot(3,1,3) plot(t,T_parameters); title('Time signal from parameters'); xlabel('T [s]'); ylabel('A');

Answer 1

[David Goodmanson]

Hi Roderick, I see that what I said about dividing your angles by a factor of 2*length(...) was incorrect because that factor was inside the argument of the angle function. So I will start over. [1] I believe your time array has one too many points. For true periodicity, the value of T_addition should not be the same at the end point as at the first point. Going with something like

N  = 1000;
dt = .01
Max_t = N*dt;
df = 1/Max_t;
t = (0:N-1)*dt
f = (0:N-1)*df

gives single-point amplitudes at the desired frequencies, zero everywhere else. The extra time point messes that up. If you try plotting the fft either way you will see a significant difference. [2] Sometimes phase variations make good sense, such as with filters, but there is not exactly such a thing as a phase spectrum in all cases. If you have tiny amplitudes between peaks, down around 1e-10 or whatever, you will still get phase angles, which are basically meaningless. [3] Point 1 in the freq array is zero frequency, point 2 matches up with the end point, 3 with end-2 etc. Your code is close but in the last for loop the the -2 terms should be +2. [4] For the fourier transform or fft of a real function, the phase angle at negative frequency is minus the phase angle at positive frequency. So in the last line of the last for loop the = phase should be = - phase. Once those changes are made the plots are reasonably good. The phases from the fft look like noise but almost all those phases belong to minuscule amplitudes; I believe the phases at the actual signal points are correct.

[Roderick]

David, thanks a lot for your detailed answer, it really helped me.

[Oleksandr Radomskyi]

Hi David, Roderick, thanks a lot for sharing this!