Impulse response from frequency response in Matlab

Question

I am trying to design equalization filter and therefore I want to define my own amplitude and phase response and then to obtain the impulse response of the filter. I thought that the output of the function freqz is the same as abs(freqz(imp)).exp(1iphase(freqz(imp))). As I want to estimate the coefficients of equalization filter I would like to get the impulse response of 1/freqz(coeff_channel), but when I try coeff_equal = real(ifft( 1./freqz(coeff_channel))) and then try finding the freqz(coeff_equal) I don't get the same results as when I define separately abs(1./freqz(coeff_channel)) and phase (1./freqz(coeff_channel)). I am attaching below the mentioned code. My question is, what is the correct way to go from the frequency response to impulse response?

clear all; close all; clc;
%%
Fs = 306.700*1e+3; % sampling rate 
Fc = 100*1e+3; % central frequency
B = 4*1e+3;  % bandwidth
N = 3072;
T = (N-1)/Fs;
t = 0:1/Fs:T; %timeccc cc

hh = [0.3 0.5 0.9 -0.1]; % Channel
%hh = [0.9 0.9 0.5 0.9]; % Channel 
figure
freqz(hh)
title('Channel')

figure
stem(hh)
title('Impulse response of channel')


baseband_signal = chirp(t,-B/2,t(end),B/2)-i*chirp(t,-B/2,t(end),B/2,'linear',90); 

% Modulation
modulated_signal = baseband_signal.*exp(i*2*pi*Fc.*t);
modulated_signal = real(modulated_signal); 
 
% Transmitted signal
tx = modulated_signal.*1.25/max(modulated_signal);
  
% Received signal 
%rx = [zeros(1,200) tx+[zeros(1,200) tx(1:(end-200))]]; 
rx = filter(hh,1,tx);
%data = xlsread('chirp100kHz_B4kHz_1plusminus.csv');
%rx = data(:,2)';

t_received = 0:1/Fs:(length(rx)-1)/Fs;

% Computing equalization filter: rx and tx swapped!
 coeff_equal = channel(tx,rx,4); % impulse response of equalization filter
 coeff_channel = channel(rx,tx,4); % impulse response of channel
 
%coeff_equal = real(ifft(abs(1./freqz(coeff_channel)).*exp(1i*(phase(1./freqz(coeff_channel))))));

% Equalization filtering
 rx_corrected = filter(coeff_equal, 1,  rx);
 
 figure
freqz(coeff_equal,1)
title('Equalization filter')

 figure
 freqz(coeff_channel,1)
title('Channel estimation')
[h,w] = freqz(coeff_channel,1);
w = w./pi;

 figure
subplot(2,1,1)
hold all;
plot(abs(freqz(coeff_channel))); 
plot(abs(freqz(hh))); 
subplot(2,1,2)
hold all; 
plot(phase(freqz(coeff_channel))); 
plot(phase(freqz(hh))); 
legend({' Channel estimation','Channel'})

 figure
subplot(2,1,1)
plot(w,20*log10(abs(1./freqz(coeff_channel)))); 
hold all;
plot(w,20*log10(abs(freqz(coeff_channel)))); 
plot(w,20*log10(abs(freqz(coeff_channel)))+20*log10(abs(1./freqz(coeff_channel)))); 
subplot(2,1,2)
hold all; 
plot(w,(phase(1./freqz(coeff_channel))).*180/pi); 
plot(w,(phase(freqz(coeff_channel))).*180/pi); 
plot(w,(phase(freqz(coeff_channel))).*180/pi+(phase(1./freqz(coeff_channel))).*180/pi); 
legend({'Equalization filter',' Channel estimation','Combined'})


 figure
subplot(2,1,1)
plot(w,20*log10(abs(freqz(coeff_equal)))); 
hold all;
plot(w,20*log10(abs(freqz(coeff_channel)))); 
plot(w,20*log10(abs(freqz(coeff_channel)))+20*log10(abs(freqz(coeff_equal)))); 
subplot(2,1,2)
hold all; 
plot(w,(phase(1./freqz(coeff_equal))).*180/pi); 
plot(w,(phase(freqz(coeff_channel))).*180/pi); 
plot(w,(phase(freqz(coeff_channel))).*180/pi+(phase(freqz(coeff_equal))).*180/pi); 
legend({'Equalization filter',' Channel estimation','Combined'})

figure
plot(t_received,rx,t_received,rx_corrected)
xlabel('Time (s)');
ylabel('Amplitude');
grid on;
legend('rx','rx_corrected')


figure
plot(t,tx,t_received,rx_corrected)
xlabel('Time (s)');
ylabel('Amplitude');
grid on;
legend('tx','rx_corrected')

figure
plot(t,tx,t_received,rx)
xlabel('Time (s)');
ylabel('Amplitude');
grid on;
legend('tx','rx')

function coeff = channel(tx,rx,ntaps)
    % Determines channel coefficients using the Wiener-Hopf equations (LMS Solution)
    % TX = Transmitted (channel input) waveform, row vector, length must be >> ntaps 
    % RX = Received (ch output) waveform, row vector, length must be >> ntaps 
    % NTAPS = Number of taps for channel coefficients
    % Dan Boschen 1/13/2020

    tx= tx(:)';   % force row vector
    rx= rx(:)';   % force row vector
    depth = min(length(rx),length(tx));
    A=convmtx(rx(1:depth).',ntaps);
    R=A'*A;       % autocorrelation matrix
    X=[tx(1:depth) zeros(1,ntaps-1)].';
    ro=A'*X;      % cross correlation vector
    coeff=(inv(R)*ro);   %solution
end

Answer 1

Equalizing a channel by inverting its frequency response is NOT recommended. This is because only systems that have a minimum phase response (all zeros inside the unit circle) have a stable causal inverse. For example, any channel with leading and trailing echos (most wireless channels) are a mixed phase system so therefore cannot simply be inverted.

The function channel that the OP used will work to determine the equalizer instead of the channel itself by swapping tx and rx (as well as adjust the relative delay between the two such that the resulting equalizer impulse response is centered within the resulting equalizer filter). This approach of that function is further detailed here as the function equalize:

Compensating Loudspeaker frequency response in an audio signal

Using tx and rx as included by the OP solves for the channel since tx convolves with the channel to give rx. By swapping them you instead get the equalizer since rx convolved with the equalizer recovers tx.

I pasted the slides below for another response but can't find that so will duplicate them here as well since they are relevant. A recursive equalizer refers to an IIR filter (feedback). Poles outside the unit circle is unstable for causal systems.