HPF doesn't work in time-domain convolution

6 visualizzazioni (ultimi 30 giorni)
Tianjia
Tianjia il 21 Feb 2023
Commentato: Tianjia il 21 Feb 2023
This is my input signal:
t=0:t_sample:1e-3;
xraw=cos(2*pi*1e6*t)+cos(2*pi*1e7*t)+cos(2*pi*1e8*t);
k=15;
N=2^k;
x=xraw(1:N); %I made my signal have a length of 2^k just to make it easier for fft and nothing else
X=fft(x)/length(x);
f=Fs/N*(-N/2:N/2-1);
and I designed a very simple High Pass Filter:
R=100;
C2=1e-13;
HPF_sym=R*C2*j*2*pi*f./(R*C2*j*2*pi*f+1);
and I want to see if multiplication in frequency domain and convolution in time-domain have the same result, so I did ifft, before I ifft I mapped the freuqency response to match with input spectrum X. ( I'm not sure if this is necessary)
for i = 1:N/2
HPF(i)=HPF_sym(i+N/2);
end
for i=N/2+1:N
HPF(i)=HPF_sym(i-N/2);
end
and then I did ifft and do convolution :
hpf=ifft(HPF);
x1=conv(x,hpf,'same');
and I want to see whether the spectrum is correct:
X1=fft(x1)/length(x1);
but X1's spectrum shows similar amplitude for pulses at 1e6 1e7 1e8 frequency, which is incorrect. Dirrect Multiplication of X and HPF_sym gives the correct spectrum. I don't know why the time-domain convolution gives the wrong answer.
I've also tried not using direct ifft to derive the time domain impulse response h(t) from H(w). Since ks/(ks+1) in s-domain corresponds to hpf(t)=dirac(t) - k * exp(-k), I've also tried using this as the time domain impulse response to convolute with x:
k=R*C2;
hpf(t) = -k*exp(-k*t);
hpf(1) = 1-k; %In order to avoid using dirac(x) function in MATLAB, I separated the h(t) function.
hpf = hpf(1:N) %Just to make sure hpf has the same length as input signal x
x2 = conv(x,hpf,'same');
X2=fft(x2)/length(x2);
semilogx(abs(X2));
And the spectrum of X2 is similar to X1, which is not the expected "high frequencies maintained and low frequencies attenuated" form. Rather, Both X1 and X2 shows a spectrum in which high frequencies and low frequencies has similar amplitudes.
I even tried the following codes, where I used ifft(X) and ifft(HPF) to do convolution in time-domain:
subplot(1,2,1);
semilogx(abs(X.*HPF));
subplot(1,2,2);
semilogx(abs(fft(conv(ifft(X),ifft(HPF),'same'))));
and the two results are still far from each other. New to signal processing and I'm really confused.

Risposte (1)

Askic V
Askic V il 21 Feb 2023
Please have a look at this code. This is to confirm that multiplication in freq. domain and convolution in time domain will provide the same result.
t = linspace(-1,1,100);
y1 = sin(t)+2*cos(2*t);
y2 = 3*sin(t)+4*cos(3*t);
L1 = numel(y1);
L2 = numel(y2); % In this case L1 is the same length as L2
% Convolution of signals
yc = conv(y1,y2);
% Now calculate via multiplication in freq. domain
L = L1+L2-1; % total length of the signal
t = linspace(-1,1,L1+L2-1); % new length of t
y1 = [y1, zeros(1,L2-1)];
y2 = [y2, zeros(1,L1-1)];
Y1 = fft(y1);
Y2 = fft(y2);
Y = Y1.*Y2;
y = ifft(Y);% calculate inverse FFT
% plot
plot(t, yc)
hold on
plot(t, y)
legend('Convolution', 'Inverse FFT');
In your case, I would calculate FFT of input signal and use not only abs (amplitude spectrum), but complete information i.e. complex numbers (they contain both amplitude and phase information), then multiply with HPF and do an inverse FFT of that.
Then compare this to convolution of signal x in time domain and impulse response of HP filter.
  3 Commenti
Askic V
Askic V il 21 Feb 2023
The problem I see is that you need an impulse response of the HPF. The real question is how to calculate impulse response which should be an array of double numbers (in time domain) and it should not have complex numbers.
Tianjia
Tianjia il 21 Feb 2023
Do you mean directly using ifft(H(w)) doesn't work? I thought about this when I found that abs(fft(conv(ifft(X),ifft(HPF))) and abs(X.*HPF) doesn't have the same result.
But I also tried using the theoretic impulse response h(t)=diract(t) - k exp(-k*t) and convolute x with this h(t) and still the result is incorrect.

Accedi per commentare.

Categorie

Scopri di più su Fourier Analysis and Filtering in Help Center e File Exchange

Prodotti


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by