Phase shift of transfer function using MATLAB

I'm trying to compute the phase shift of a transfer function (accelerance) in MATLAB, but I'm having issues.

The transfer function was computed by dividing the magnitude of the system output by the input (please see here).

I've been using pwelch() to compute all of the functions in the frequency domain, but from the best of my understanding, the Welch PSD estimate outputs a real-valued vector, so there's no way to get a phase shift.

When I try to use the fft() function, the result (magnitude) doesn't look right, so I'm not sure if I would trust the phase shift. Below is the code I used to compute the transfer function:

% Welch Acceleration
[pa,f] = pwelch(ay,1000,10,0:10:10000,Fs);
% Welch Force
[pf,f] = pwelch(fy,1000,10,0:10:10000,Fs);
% Transfer function
H = pa./pf;


If someone could point me in the right direction it would be appreciated.
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.6548966765403748, "perplexity": 509.5753185490501}, "config": {"markdown_headings": false, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2021-17/segments/1618039610090.97/warc/CC-MAIN-20210422130245-20210422160245-00099.warc.gz"}