Need help to solve Complex Demodulation problem

14 views (last 30 days)
I am trying to do Complex Demodulation. I have noisy signal y.
My task is extract amplitude and phase at frequency of 20 Hz and 30 Hz from y.
In this plot you can see my result of amplitude demodulation at 20 Hz.
My problem is -- after low pass filtering amplitude of the demodulated signal (19) becomes 2 times higher than original signal y (40).
I want the amplitude of demodulated signal to be 19 like in the original data y.
There is something wrong with low pass filter or modulation procedure.
Here is code:
clc
clear
fs = 100; % sampling frequency (Hz)
Ts = 1/fs;
f1 = 10; % Hz
f2 = 20;
f3 = 30;
t=(1:1000)*Ts; % time
t=t';
noise = (30*rand(1000, 1));
y = 10.*sin(2*pi*f1*t)+20*sin(2*pi*f2*t)+30*sin(2*pi*f3*t) + noise; % Noisy input signal
s = abs(fft(y)./500);
f = (1:length(y))/length(y).*fs;
f=f';
mod=y.*(exp(-1i*2*pi*f2*t)); % complex demodulation
[b,a] = butter(10,f2/(fs/2),'low'); % Butterworth filer
s4a = filtfilt(b,a,imag(mod));
s4b = filtfilt(b,a,real(mod));
s4 = sqrt(s4a.^2+s4b.^2); % Demodulated Amplitude
phase = atan2(s4a, s4b); % Demodulated Phase
subplot(3,1,1);
plot(t,y)
title('Noisy input signal');
xlabel('t')
subplot(3,1,2);
plot(f,s)
title('FFT Spectrum of Input signal');
xlabel('f, Hz')
subplot(3,1,3);
plot(t,s4)
title('Amplitude Demodulated signal at 20 Hz');
xlabel('t')

Answers (1)

Star Strider
Star Strider on 17 Nov 2021
My problem is -- after low pass filtering amplitude of the demodulated signal (19) becomes 2 times higher than original signal y (40).
It really is not.
Here, I changed the first subplot so that ylim went from 0 to the previous maximum value of ylim, and that demonstrates that the amplitude is actually times higher than the demodulated signal, so the demodulated signal amplitude is actually 66% of the original signal. That can be seen in the RMS calculations (added at the end of the code).
This is of course to be expected, since much of the signal energy was removed in the filtering process.
fs = 100; % sampling frequency (Hz)
Ts = 1/fs;
f1 = 10; % Hz
f2 = 20;
f3 = 30;
t=(1:1000)*Ts; % time
t=t';
noise = (30*rand(1000, 1));
y = 10.*sin(2*pi*f1*t)+20*sin(2*pi*f2*t)+30*sin(2*pi*f3*t) + noise; % Noisy input signal
s = abs(fft(y)./500);
f = (1:length(y))/length(y).*fs;
f=f';
mod=y.*(exp(-1i*2*pi*f2*t)); % complex demodulation
[b,a] = butter(10,f2/(fs/2),'low'); % Butterworth filer
s4a = filtfilt(b,a,imag(mod));
s4b = filtfilt(b,a,real(mod));
s4 = sqrt(s4a.^2+s4b.^2); % Demodulated Amplitude
phase = atan2(s4a, s4b); % Demodulated Phase
subplot(3,1,1);
plot(t,y)
title('Noisy input signal');
xlabel('t')
grid
ylim([0 max(ylim)])
subplot(3,1,2);
plot(f,s)
title('FFT Spectrum of Input signal');
xlabel('f, Hz')
grid
subplot(3,1,3);
plot(t,s4)
title('Amplitude Demodulated signal at 20 Hz');
xlabel('t')
grid
Original_RMS = rms(y)
Original_RMS = 31.4731
Filtered_RMS = rms(s4)
Filtered_RMS = 20.8819
Filtered_RMS/Original_RMS
ans = 0.6635
.
  4 Comments
Star Strider
Star Strider on 18 Nov 2021
Sure!
For a filter with a passband of 20 to 30 Hz, this works —
Fs = 500; % Use Correct Sampling Frequency
Fn = Fs/2;
Wp = [20 30]/Fn; % Passband Frequency (Normalised)
Ws = [15 35]/Fn; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 60; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
figure
freqz(sos, 2^16, Fs) % Filter Bode Plot
% signal_filt = filtfilt(sos, g, signal); % Filter Signal
If the ‘Wp’ and ‘Ws’ are vectors, the MATLAB filter functions define a bandpass filter by definition, and if the last argument to the ellip call is 'stop' will design as bandstop filter instead. (I left the filtfilt call commented-out so it wouldn’t throw an error when I ran the code here.)
Experiment with the stopbands to get different results.
.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!