How to calculate magnitude/phase by using a s-function?

Question

Kevin 2012-12-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/56524-how-to-calculate-magnitude-phase-by-using-a-s-function

Hey guys,

I'm trying to calculate the magnitude and the phase of complex input signals by using a s-function (level-1 m-file). There are 3 different input signals in the form z = a+b*j. Both, the real and the imaginary parts of each signal, are merged to a bus signal. So the s-function has 6 inputs in total (3x real + 3x imag). Now, I'd like to determine the magnitude and phase of the corresponding pairs of inputs and put them on the bus again (3x magnitude + 3x phase). However, when compared to the "Real-Imag to Complex2 and "Complex to Magnitude-Angle" blocks in Simulink there are differences of up to 6 rad considering the phase!! Why's that? Do you have any advice of what might went wrong?

To be able to follow my thoughts, the s-function source code is attached:

function [sys,x0,str,ts,simStateCompliance] = test_reimg(t,x,u,flag)
switch flag,
%%%%%%%%%%%%%%%%%%
% Initialization %
%%%%%%%%%%%%%%%%%%
  case 0,                                                
    [sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes;    
%%%%%%%%%% 
% Update %
%%%%%%%%%%
  case 2,                                               
    sys = mdlUpdate(t,x,u);
%%%%%%%%%%
% Output %
%%%%%%%%%%
  case 3,                                               
    sys = mdlOutputs(t,x,u);    
%%%%%%%%%%%%%
% Terminate %
%%%%%%%%%%%%%
  case 9,                                               
    sys = [];
    otherwise
      DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
  end
%end sfundsc2
%
%    =============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
global n;
cs = 3;              % Anzahl der Eingangssignale der Form y = a+bj
n = cs*2;           % Gesamtanzahl der Zustände für Realteil und Imaginärteil
sizes = simsizes;
sizes.NumContStates  = 0;
sizes.NumDiscStates  = n;
sizes.NumOutputs     = n;
sizes.NumInputs      = n;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0  = [0 0 0 0 0 0];
str = [];
ts  = [-1 0;]; % Sample period 
% speicfy that the simState for this s-function is same as the default
simStateCompliance = 'DefaultSimState';
% end mdlInitializeSizes
%
%=======================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=======================================================================
%
function sys = mdlUpdate(t,x,u)
sys = x;
%end mdlUpdate
%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u)
global n;
r = [];
phi = [];
limit = n/2
for x = 1:limit
  r(x) = abs(u(x)+u(x+limit)*1i)
  phi(x) = angle(u(x)+u(x+limit)*1i)-pi/4
    %{ 
    r(x) = sqrt(u(x)^2+u(x+limit)^2);
     if u(x) > 0    
        phi(x) = atan(u(x+limit)/u(x))  - pi/4              % 1./4. Quadrant
    elseif u(x) < 0  && u(x+limit) > 0                       % 2. Quadrant
        phi(x) = atan(u(x+limit)/u(x)) + pi - pi/4
    elseif u(x) < 0  && u(x+limit) < 0                       %3. Quadrant
        phi(x) = atan(u(x+limit)/u(x)) - pi - pi/4  
    end
   %} 
  end
sys = [r phi];
%end mdlOutputs