find function not working with double inequality

Question

S 2024-4-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2105771-find-function-not-working-with-double-inequality

编辑： S 2024-4-21

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Voss 2024-4-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2105771-find-function-not-working-with-double-inequality#answer_1439786

编辑：Voss 2024-4-11

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 )

10 个评论
显示 8更早的评论隐藏 8更早的评论

S 2024-4-15

编辑：S 2024-4-15

@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 2024-4-15

在 MATLAB Online 中打开

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

Voss 2024-4-18

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 2024-4-21

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

请先登录，再进行评论。

find function not working with double inequality

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

find function not working with double inequality

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论