find function not working with double inequality

I am trying to use a double inequality to define in my loop which outputs to display together from my filter. I am not getting any error, but my figure 2 is showing all outputs instead of just the few which fall into the designated range. To check this specific case I looked into my array and found that figure 2 should be only displaying only output 5, 6, and 7. I am redefining my output variable to include indf which is supposed to find those outputs as they fall in that range. I am unsure what is going wrong. Thank you for your time!

close all
clear all
clc
%
fs = 20e3; 
numFilts = 32; % 
filter_number = numFilts;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
% input signal definition can increase Nperiods instead of zero padding for
% resolution
signal_freq = 1000; % Hz
Nperiods = 15; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
t = linspace(0,Nperiods/signal_freq,200*Nperiods); % 
input = sin(2*pi*signal_freq*t) + 0*rand(size(t));
figure%2
hold on
for ii = 1:filter_number
    IR = gammatone(order, CenterFreqs(ii), fs);
    [tmp,~] = freqz(IR,1,1024*2,fs);
    % scale the IR amplitude so that the max modulus is 0 dB
    a = 1/max(abs(tmp));
    IR_array{ii} = IR*a; % scale IR and store in cell array afterwards
    [h{ii},f] = freqz(IR_array{ii},1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(IR_array{ii})
end
title('Impulse Response');
hold off
xlim([0 1000]);
%% 3D
figure
hold on
b=50; %bandwidth
indf = find(signal_freq-b<CenterFreqs<signal_freq+b); 
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);
    
    plot3(t,ii*ones(size(t)),output(ii,:)) 
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
view(-40,60); % view angles 
% display integer ticks on Y axis
ax = gca;
ax.YTick = unique(round(ax.YTick));
output_sum = sum(output,1);
% plot(t,output_sum)
% LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals');
xlabel('Time (s)');
ylabel('Output #');
zlabel('Amplitude');
hold off

Answer 1

[Voss]

This is not the correct way to check multiple inequalities:

indf = find(signal_freq-b<CenterFreqs<signal_freq+b);

Instead do:

indf = find( signal_freq-b<CenterFreqs & CenterFreqs<signal_freq+b );

or:

indf = find( abs(CenterFreqs-signal_freq) < b )

[Voss]

find(signal_freq-b<CenterFreqs<signal_freq+b)

is interpreted as

find( (signal_freq-b<CenterFreqs) < signal_freq+b )

so

signal_freq-b<CenterFreqs

is done first, and then the result of that, which is a logical array (i.e., an array of true and/or false values), is compared to signal_freq+b. It's the same as if you said:

temp = signal_freq-b<CenterFreqs;

find( temp < signal_freq+b )

When comparing a logical (true or false) value to a numeric value, true is considered 1 and false is considered 0, so you're comparing ones and/or zeros to signal_freq+b, which is not what you meant to do.

[Steven Lord]

FYI if you edit the user's code in MATLAB Editor, Code Analyzer flags the line where that comparison occurs with a shorter version of the argument you made in your comment as explanation. I don't remember off the top of my head when that Code Analyzer warning was introduced, but I believe it's been around for several years now.

[S]

@Voss oh! That makes sense! Thank you!

So I tried adding another loop right after to make a second 3d plot:

figure

hold on

s=200;

indf = find( signal_freq-s<signal_freq-b & signal_freq-s<signal_freq+b );

for ii = 1:numel(indf)

output(ii,:) = filter(IR_array{indf(ii)},1,input);

plot3(t,ii*ones(size(t)),output(ii,:))

LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each

end

view(-40,60); % view angles

% display integer ticks on Y axis

ax = gca;

ax.YTick = unique(round(ax.YTick));

output_sum_to_subtract = sum(output,1);

new=output_sum-output_sum_to_subtract;

plot3(t,ii*ones(size(t)),new)

LEGs{ii+1} = ['Sum'];

legend(LEGs{:})

legend('Show')

title('Filter output signals');

xlabel('Time (s)');

ylabel('Output #');

zlabel('Amplitude');

hold off

But I am getting:

Incorrect dimensions for matrix multiplication. Check that the number of

columns in the first matrix matches the number of rows in the second matrix.

To operate on each element of the matrix individually, use TIMES (.*) for

elementwise multiplication.

Error in t4112 (line 104)

plot3(t,ii*ones(size(t)),new)

But when I check both new and t are 1x3000?

[Voss]

You're welcome!

This doesn't seem right:

indf = find( signal_freq-s<signal_freq-b & signal_freq-s<signal_freq+b );

Do you think you meant to do this?

indf = find( signal_freq-s>signal_freq-b & signal_freq-s<signal_freq+b );

Anyway, each term has a signal_freq which can be subtracted out from all, leaving

indf = find( -s>-b & -s<+b );

which is the same as

indf = find( s<b & s>-b );

Is that the condition you meant?

[Voss]

Anyway, that error happens because ii is empty, which happens because indf is empty.

[S]

So I am trying to get for ex. if signal freq is 1000 and s is 100 and b is 100. With the first 3d plot I am trying to find the 900-1100 range and in this second plot I am trying to find 800 to 900 and 1100 to 1200 valued. @Voss

[Voss]

In that case, you can use something like

s=150;

indf = find( abs(CenterFreqs-(signal_freq+s)) < b | abs(CenterFreqs-(signal_freq-s)) < b );

signal_freq is 1000, so signal_freq+s is 1150 and signal_freq-s is 850, so the condition above finds CenterFreqs that are within b (50) of either of those, i.e., within 50 of 850 (which is 800 to 900) or within 50 of 1150 (which is 1100 to 1200).

close all
clear all
clc
%
fs = 20e3; 
numFilts = 32; % 
filter_number = numFilts;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
% input signal definition can increase Nperiods instead of zero padding for
% resolution
signal_freq = 1000; % Hz
Nperiods = 15; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
t = linspace(0,Nperiods/signal_freq,200*Nperiods); % 
input = sin(2*pi*signal_freq*t) + 0*rand(size(t));
figure%2
hold on
for ii = 1:filter_number
    IR = gammatone(order, CenterFreqs(ii), fs);
    [tmp,~] = freqz(IR,1,1024*2,fs);
    % scale the IR amplitude so that the max modulus is 0 dB
    a = 1/max(abs(tmp));
    IR_array{ii} = IR*a; % scale IR and store in cell array afterwards
    [h{ii},f] = freqz(IR_array{ii},1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(IR_array{ii})
end
title('Impulse Response');
hold off
xlim([0 1000]);

%% 3D
figure
hold on
b=50; %bandwidth
% indf = find(signal_freq-b<CenterFreqs<signal_freq+b); 
indf = find( abs(CenterFreqs-signal_freq) < b );
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);
    plot3(t,ii*ones(size(t)),output(ii,:)) 
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
view(-40,60); % view angles 
% display integer ticks on Y axis
ax = gca;
ax.YTick = unique(round(ax.YTick));
output_sum = sum(output,1);
% plot(t,output_sum)
% LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals');
xlabel('Time (s)');
ylabel('Output #');
zlabel('Amplitude');
hold off

%% perform FFT of signal : 
[freq,fft_spectrum] = do_fft(t,output_sum);
figure
plot(freq,fft_spectrum,'-*')
xlim([0 1000]);
title('FFT of output sum signal')
ylabel('Amplitude');
xlabel('Frequency [Hz]')

% second 3d plot
figure
hold on
s=150;
indf = find( abs(CenterFreqs-(signal_freq+s)) < b | abs(CenterFreqs-(signal_freq-s)) < b );
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);         
    plot3(t,ii*ones(size(t)),output(ii,:)) 
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
view(-40,60); % view angles 
% display integer ticks on Y axis
ax = gca;
ax.YTick = unique(round(ax.YTick));
output_sum_to_subtract = sum(output,1);
new=output_sum-output_sum_to_subtract;
plot3(t,ii*ones(size(t)),new) 
LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals');
xlabel('Time (s)');
ylabel('Output #');
zlabel('Amplitude');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)
time = time(:);
data = data(:);
dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length
%% use windowing or not 
% no window
fft_spectrum = abs(fft(data))*2/nfft;
%fft_spectrum = abs(fft(data.*window),nfft))/sum(window);
% % hanning window
window = hanning(nfft);
window = window(:);
fft_spectrum = abs(fft(data.*window))*4/nfft;
% one sidded fft spectrum  % Select first half 
if rem(nfft,2)    % nfft odd
    select = (1:(nfft+1)/2)';
else
    select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IR = gammatone(order, centerFrequency, fs)
% Design a gammatone filter
earQ = 9.26449;
minBW = 24.7;
erb = ((centerFrequency/earQ)^order + (minBW)^order)^(1/order);% we use the generalized
% your code
B = 1.019 * erb;
f = centerFrequency;
tau = 1. / (2. .* pi .* B);
tmax = tau + 20/(2*pi*B);
t = (0:1/fs:tmax);
temp = (t - tau) > 0;
IR = ((t - tau).^(order - 1)) .* exp(-2*pi*B*(t - tau)).* cos(2*pi*f*(t - tau)) .* temp;
end

[S]

@Voss This makes so much sense! Thank you! So I'm just trying to better understand the inequality and how to write it correctly.

If I did something like:

indf = find(abs(CenterFreqs-(signal_freq-s)) < b );

Instead in our new segment would this now just be finding the filter with the 150 lower value?

While if I included the other half of the find it would be finding the filter within 150 higher right?

Like:

indf = find( abs(CenterFreqs-(signal_freq+s)) < b );

[Voss]

indf = find(abs(CenterFreqs-(signal_freq-s)) < b );

finds elements of CenterFreqs that are within b (=50) of signal_freq-s (=850), so within the range 800 to 900, exclusive.

indf = find( abs(CenterFreqs-(signal_freq+s)) < b );

does the same within the range 1100 to 1200, exclusive.

So what you said is right, except it's 100 (which is 2*b) not 150 (which is s).

[S]

@Voss oh makes sense! So the first line is for lower freqs, while the second is for higher freqs. Thank you!