Need help to solve Complex Demodulation problem

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')

Answer 1

[Star Strider]

‘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 1½ 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

.

[Steve Ts]

@Star Strider Thank you for your answer. But I'm still unclear how I can achieve the result I want.

I need to extract this part of the signal - 20*sin(2*pi*f2*t) (at a frequency of 20 Hz with an amplitude of 20) from the input signal y.

[Star Strider]

To extract one frequency or a band of frequencies, use the bandpass function (introduced in R2018a) and choose 'ImpulseResponse','iir' to have it design an effficient elliptic filter with the chosen passband. To extract the 20 Hz signal, set the passband to [15 25] for best resullts. Passbands that are too narrow can produce an unstable filter, regardless of the filter design, and there is always spectral leakage due to sampling and other effects.

If only earlier MATLAB releases are available, designing an equivalent filter is straightforward. I can help with that if necessary.

.

[Steve Ts]

@Star Strider could i ask you to help me with complex demodulation extraction using Passband-filter? I really need it.

[Star Strider]

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.

.