texts: Number of samples in one period of subcarrier after IFFT

texts:
I have a quick question about the number of samples in one period of a sub carrier in OFDM systems. The IFFT process generates N-parallel sub carriers in discrete time domain which are sent simultaneously over the sub carrier frequencies.

What is the number of samples in one period of a sub carrier after IFFT?

Thanks for your reply. Let me give you more details about my query with an illustration as below. The picture shows the sub carriers as analog but lets take them as discrete before DAC and CP addition. The system has 4 sub carriers and hence they have different frequencies. My questions is here at this point; again assuming they are just the output of the IFFT and digital-
Does each sub carrier has 4 samples in every cycle? So lets say the duration of OFDM symbol is T seconds. The subcarrier1 has 1 cycle and I am visualizing that the sub carrier has 4 samples during the T seconds.
Subcarrier2 has 2 cycles (period) and hence it has 2x4=8 samples during the OFDM T seonds.
And for subcarrier3 it has 3x4=12 samples.
For subcarrier4, there are 4x4=16 samples during the OFDM symbol (T seconds length).

And in total, the OFDM symbol has 36 samples in the symbol length T seconds.

Am I correct interpreting this?


texts: Design FIR filter with half zero coefficions

texts:
I noticed a lot of hardware IC use special FIR filters with half zero of coefficients to minimize multiplication operations, especially in decimators/interpolators. E.g.
[ 6 0 -19 0 47 0 -100 0 192 0 -342 0 572 0 -914 0 1409 0 -2119 0 3152 0 -4729 0 7420 0 -13334 0 41527 65536 41527 0 -13334 0 7420  0 -4729 0 3152 0 -2119 0 1409 0 -914 0 572 0 -342 0 192 0 -100 0 47 0 -19 0 6]
or
[-12 0 84 0 -336 0 1006 0 -2691 0 10141 16384 10141 0 -2691 0 1006 0 -336 0 84 0 -12]
I took these coefficients from IC datasheet. I see the center coefficient is 0.5 and other odd coefficients are zero. It effectively halves the mul operations.
But I can't find what the type of the FIR filter is. Is there any good paper? How to create a filter with desired number of taps with this feature?


texts:
This is a difference equation to a causal LTI system:
$y[n] = ay[n - 1] + x[n] - a^Nx[n - N]$
Where N is a positive integer. I need to determine the impulse response of the system, so I have the equation:
$h[n] = ah[n-1] + \delta[n] - a^N\delta[n - N]$
h[n] is simple to find for the case where n < 0, n = 0, or n = 1. However, the response could be different for the case n = 2, depending on the value of N (in this case, N = 1 or N = 2):
$h[2] = ah[1] + 0 - a^1\delta[2-1] = a^2$
$h[2] = ah[1] + 0 - a^2\delta[2-2] = 0$
So which one is it? Intuition tells me the top one as the impulse response will be decaying, but we only know that N is positive.


texts:
I have an image depicting a results of gel electrophoresis (it does not really matter - it just means: no bodies, faces, predictible objects) and I have a strong suspicion that a part of the image background was copied (probably Photoshop Clone Stamp) to cover an unwanted peak. The background is really just noise, but it should be almost perfect copy. I am not entirely sure how big is the copied area so just straightforward cross-correlation with an image part won't work.

How do I prove that? I prefer Matlab for such tasks.


texts: How to determine whether a speech segment is voiced/unvoiced?

texts:
I want to determine whether a speech frame is voiced/unvoiced. Out of many methods found while searching, one method said find energy of the frame and if it is above a certain threshold, mark it as voiced. Now, my question is how should I determine this 'threshold value'? Is it by trial and error or are there any set of rules?

In my attempts, I resorted to a simple idea of looking at the energy plots and setting a threshold value accordingly. It served me well, but I want to know whether it was just beginner's luck.


texts: Equation visualization of complex sinusoid function

texts:
I'm trying to learn Fourier Transform & Signal

I can visualize how this expression $ A \sin(2 \pi f t)$ could turn into that curvy sinusoidal signal.
But i couldn't visualize how $F(\nu)e^{2 \pi i f t}$ would look like as a signal....

Why the natural exponent suddenly popped up? what does it do?
i understand that $F(\nu)$ gives information of the amplitude and the phase shift, but what sort of information i could get from $e$ ?
Does all signals always decaying/growing by time?

I saw this a lot in Laplace Transform: $e^{-st}$ but don't understand how it would look like as a signal.


texts: Why do some radio stations still use amplitude instead of frequency modulation?

texts: Why is frequency modulation (FM) not used on the current amplitude modulation (AM) radio stations? The Police have FM 2 way radios that use only 5 kilocycles of modulation. Instead of the current AM radio stations that use the same bandwidth of 5 kilocycles. The limiter & discriminator or ratio detector filter out the static & loud pops that are on AM; Thus the FM would have much less static & noise than AM. FM channels would occupy the same bandwidth as AM channels. For example if WABC in New York used frequency modulation over the current frequency of 770 they would need to have the oscillator changed but the power amp would stay the same. But the receiving radios would need to be replaced with a radio to demodulate frequency modulation at 770.


texts:
i'm trying to implement an arbitrary IIR filter design. I'm decomposing all the algorithm into smaller pieces that i can test.

I currently am on the part where it designs a FIR filter from an arbitrary response.

Matlab's main function is here :

function h = lslevin (N ,om, D, W)
% Complex Least Squares FIR filter design using Levinson' s algorithm
%
% h filter impulse response
% N filter length
% om frequency grid (0 <= om <= pi)
% D complex desired frequenccy response on the grid om
% W positive weighting function on the grid om
om = om(:);
D = D(:);
W = W(:);
L = length(om);
DR = real(D) ;
DI = imag(D);
 a = zeros(N,1);
 b = a;
% Set up vectors for quadratic objective function
% (avoid building matrices)
dvec = D;
evec = ones(L, 1);
e1 = exp(j*om);
for i=1:N,
    a(i) = W.'*real (evec);
    b(i) = W.'*real (dvec);
    evec = evec.*e1;
    dvec = dvec.*e1;
end
a=a/L;
b=b/L;
% Compute weighted l2 solution
h = levin(a, b) ;

and levin function :

function x = levin(a ,b)
% function x = levin(a ,b)
% solves system of complex linear equations toeplitz (a)*x = b
% using Levinson's algorithm
% a ... first row of positive definite Hermitian Toeplitz matrix
% b ... right hand side vector
%
% Author: Mathias C. Lang, Vienna University of Technology, AUSTRIA
% 9-97
a = a(:);
b = b(:);
n = length(a);
t = 1;
alpha = a(1);
x = b(1) / a(1);
for i = 1:n-1,
    k = -(a(i+1 :-1: 2)' * t)/alpha;
    t = [t;0] + k*flipud( [conj(t);0]);
    alpha = alpha*( 1 - abs(k)^2);
    k = (b(i+1) - a(i+1:-1: 2)' * x) /alpha;
    x = [x;0] + k*flipud(conj(t));
end

an example function can be done here (it should be a simple bandpass):

om1=pi*[linspace(0,0.23,230),linspace(0.3,0.5,200),linspace(0.57,1,430)];
D1=[zeros(1,230), exp(-j*om1(231:430) *2), zeros(1,430)]; % zeros are bandstop others bandpass
W1=[10* ones(1,230), ones(1,200), 10 *ones(1,430) ];%weight function
h = lslevin(61,om1, D1, W1) ;
wavwrite(h, 44100,  'e:\tmp2.wav'); % write the impulse response

but when i compute the Impulse Response, i have something like that :

which is far from the original design

can you help me pinpoint the errors in the code ?

thanks a lot

Jeff

EDIT1 : here is the final FFT of the IR from the answer given below:


texts: how to stop matlab from rounding my values and keep them at same precision during calculations?

texts:
here are my IIR coefficients:

CoefIIR=[1.630078934734 1.267893338385 1.112390120011]

and here are the result given to me by matlab commands:

CoefIIR =

1.6301    1.2679    1.1124

How can I keep the precision for data storage and during calculations?


texts: How to calculate Integral Non Linearity (INL) the direct way

texts:
In the situation bellow we have a 3bit Analog to Digital Converter (ADC).

where avarage the LSB from the non ideal system is

I want to calculate the Integral Non Linearity (INL) purely from the ideal and actual columns of the table shown bellow.
I have tried to calculate INL for Logic 001 ,V_ideal=0.05 V_real=0.091 avarage step=0.08
INL(001)=(0.07-0.05)/0.091=-0.071
INL(010)=(0.15-0.15)/0.091=0
I get a totally different numvers then the INL in the solution.
Where did i go wrong? i did everything like the formal straight forward definition.
Thanks.


texts: What is causing these spikes in the output of the low pass filter?

texts:
I am low pass filtering the PCM of an arbitrary Bell 103 (300 bps) signal. However, the output of the low pass filter shows random spikes in amplitude.

I have tried to increase the bandwidth, but no luck.
Can someone explain why these spikes arise? My goal is to do some sort of envelope detection on this signal afterwards, to decode the 1's and 0's. For example by using a rectifier first.

The filter coefficients:

coeffs = [0.265842424474100 0.243276902464240 0.271096721997119 0.243276902464240 0.265842424474100]

The Matlab script used:

close all;
clear all;

% Generate a random bit sequence as input
bitsIn = randi([0 1], 16, 1)';

% Generate the amplitude of the FSK signal
amplitude = randi([7000 10000], 1, 1);

% Generate a randomized length of randomized startup noise
% x = randi([-20 20], 1, randi([0 3000], 1, 1));
x = []; % TODO

% FSK module the input bit stream
for i = 1:length(bitsIn)
    if (bitsIn(i) == 1)
        freq = 1270;
    else
        freq = 1070;
    end

    % Append the new PCM data for the new bit, with random phase
    phase = (2 * rand(1)) * pi;
    x = [x (amplitude * sin((2 * pi * freq * (0:26) / 8000) + phase))];
end

% Lowpass filter the resulting FSK samples
d = fdesign.lowpass('N,Fp,Fst,Ast',4,1070,1270,10,8000);
Hd = design(d,'equiripple');
y = filter(Hd,x);

figure();
plot(x);
title('Input data');
figure();
plot(y);
title('Low pass output');
```


texts: Gradient of Total Variation of Magnitude of Complex Function for Denoising

texts:
Say I have a complex function $f^*$ (e.g. a MRI image) that has a near piece-wise constant magnitude, but a non constant phase.

If I have an optimization problem to find $f^*$ and set up an objective function with a total variation term (e.g. for denoising or compressed sensing) it usually has the following form:

$$
obj_1(f) = \ldots + \text{TV}(f)
$$

However, since I assume that $f$ has a piece-wise constant magnitude, I think it might be better to use:

$$
  obj_2(f) = \ldots + \text{TV}(|f|)
$$

However, for a gradient based solver, one would have to know the gradient of obj2. The gradient for $obj_1(f)$ is: $\text{TV}'\left(TV(f)\right)$. What is the gradient of $obj_2(f)$?

Update:

Intuitively I would assume something like the following (since the phase has no influence on $obj_2$, leave the phase "untouched"):

$$
\text{TV}'\left(TV(|f|)\right)* e^{i \arg(f)}
$$

However, my knowledge in complex analysis is very limited and I am not sure if this makes sense.


texts: For SNR calculation, should I use FFT|X+iY| or FFT|R|?

texts:
I use a demodulator that outputs X and Y quadratures, with the R and Phase which can be calculated.

For SNR measurements, when plotting the FFT to find the SNR a signal at a certain frequency, should I be plotting the FFT of |X+iY| or the FFT of R?

My understanding would be to do the FFT of R so that the phase information can be completely canceled out, and R being the amplitude (the signal strength) is what matters to the SNR calculations.

Thank you


texts: FFT Pitch Detection methods: Autocorrelation or other?

texts:
I've seen many posts here that state that an FFT is almost all you need to get the fundamental frequency from an input stream.  The problem is that, with sound, this doesn't work; we need to look for periodicity as well.  I know that autocorrelation is useful in this arena, but even then, I am unfamiliar with how one would extract the fundamental frequency from the data after autocorrelation.  Also, I am aware that there are alternative methods as well.

I am not interested in the easiest solution; I am interested in 1: the most efficient solution and 2: the most accurate solution, especially if these are mutually exclusive.


texts:
I have I and Q samples of a time domain signal, sampled at a rate T, both of length L.

I create the complex envelope, C = sampleI+j*sampleQ.

Then take the FFT of C based on Matlab's FFT help algorithm:

Fs=1/T;
NFFT = 2^nextpow2(L)
Y = fft(C,NFFT)/L
fvar = (0:NFFT/2)*(Fs/NFFT)
plot(fvar,2*abs(Y))

The frequency content of the complex envelope should only contain positive frequencies but the FFT returns results symmetric in frequency so is there anything I need to do to C or to Y to interpret the results properly?  Do I just throw away the upper half of the frequency components in the FFT result?


texts:
I'm currently working on rectifying a respiratory noisy signal shown below:

I've already tried to filter the noise as you can see in the image below (full image):

The red one is the noisy signal whereas the blue one is the one got after applying the non-linear detrending:

dt_ecgl = detrend(ecgl);
opol = 6;
[p,s,mu] = polyfit(t,ecgnl,opol);
f_y = polyval(p,t,[],mu);
dt_ecgnl = ecgnl - f_y;

And the green one is the reference signal (correct one) which I'm aiming to have as result of filtering. So my question is how can I acheive a better result? My goal is to have the blue and the green ones conbined (same exact shape) in other words how can I denoise such deformed signal in MATLAB or using any other platform?.


texts:
Can somebody please tell me why the system $y[n] = x(2n)$ is non-causal ?

I know that causal systems depend on the past and present values of input and
this system satisfies the condition. So why is the answer that the system is non-causal ?

Thank you.


texts: machine learning on multiple signal at once

texts:
I am very new signal processing and I am asking this with the expectation that someone gives me a concrete guideline on how to approach my problem.

The problem I need to solve is that I have 2 signals (shown in green and yellow) which look like this as below in original scales

and log scale

Based on the interaction (e.g. cutting each other) of these 2 signals, few features are calculated like i.e. range where both signals are low, the range where signals are up and only a threshold difference in amplitude value but now I need some kind of Machine learning may be RNN or LSTM or 1D CNN. I have a set of pre-calculated values and corresponding signals. Now I am wondering how to apply machine learning techniques while looking at both these signals together.


texts: Understanding MATLAB/Mathematica's bilinear interpolation

texts:
I am trying to get a working understanding of how to resize images using bilinear and bicubic transformations. I coded a resizing function based on things I've read here and Wikipedia, etc., but I get different results than what MATLAB and Mathematica give.

A minimal working example is the image im=[100 50 0] in MATLAB or im={1, .5, 0} in Mathematica. My algorithm (and "common sense") tells me that resizing this image to 5 pixels should give
im=[100 75 50 25 0] and i={1, .75, .5, .25, 0}, respectively. However, I get
im=[100 80 50 20 0] and i={1, .8, .5, .2, 0} from MATLAB and Mathematica when forcing the bilinear method. Specifically, the MATLAB commands

im=[100 50 0];
imresize(im,[1,5],'bilinear')

returns

ans =  100.0000   80.0000   50.0000   20.0000         0

They are obviously doing the same, different thing than my code (and common sense).

Can someone explain their algorithm to me? I am guessing it my mistakes relate to how to "index" the new image. The original image, with only one row, has a constant "y" value of 1, with columns x=1,2,3. I am assuming that the new image could be viewed as having columns indexed as x=1, 1.5, 2, 2.5, 3, and interpolation is needed to compute the values of f(1.5,1) and f(2.5,1).

Thanks.


texts: Identifying whether or not a cyclostationary signal is noisy or not using the cyclic autocorrelation

texts:
I am trying to determine whether or not a given signal has been corrupted by Gaussian noise, either bandlimited (with a filter) or not. The signal in question is a BPSK or PAM signal that is upsampled and shaped by a Nyquist filter and then frequency modulated. Due to the signal being cyclostationary, I figured the best way was to use the Cyclic Autocorrelation but I am having a hard time understanding how to code it up in Matlab. This paper shows shifting the FT by some frequency $\alpha$ and then doing $F^{-1}(X(f-\alpha)*X'(f + \alpha))$ where $F^{-1}(\cdot)$ denotes the inverse DFT operation. Would I be able to just multiply it in time by a complex exponential at frequency $\alpha$ or is there more to it? Is this a valid approach or am I on a goose chase?

Thanks


texts: Matlab: How to apply constant modulus algorithm in wireless communication

texts:
I am trying to apply the Constant Modulus Algorithm which is a blind equalization algorithm in CDMA communication. I am following the code for CMA and a great explanation given in link MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm

For different signal-to-noise ratio, I want to see the performance of CMA. But I cannot understand where to plug in the CMA part in CDMA. CMA should be at the receiver, but it is not clear to  me exactly how to do so and if I need to apply the same FIR channel that is given in the CMA code. It will be very helpful is somebody can inform / show how to apply CMA. Thank you.
The code for CDMA is

clc;
clear all;
code_length = input('Enter Code Length = ');
data_stream = input('Enter Data Bit Stream = ');
%Generating Walsh Code
code = [-1 -1; -1 +1];
[r1 c1] = size(code);
disp('Walsh code generated = ');
while r1<code_length
    code = [code,code;code,-1*code];
    [r1 c1] = size(code);
    disp(code);
end
code_length = length(code); %length (number of bits) of code
Y = size(data_stream);
N = Y(1); %number of unique senders / bit streams
I = Y(2); %number of bits per stream
T = []; %sum of all transmitted and encoded data on channel
data_received = []; %vector of reconstructed bits at receiver
%show data bits and codes
disp('Data bit stream to be transmitted = '),disp(data_stream);
disp('Walsh codes used for transmission = '),disp(code);
%encode bits and transmit
G = zeros(I,code_length);
for n = 1:N
    Z = zeros(I,code_length);
    for i = 1:I
        for m = 1:code_length
            Z(i,m) = [data_stream(n,i)*code(n,m)];
        end
    end
    G = G + Z;
end
%show channel traffic
for i = 1:I
    T = [ T G(i,:) ];
end
disp('Resulting traffic on the channel = '),disp(T);
%decode and reconstruct
for n = 1:N
    tot = zeros(1,I);
    R   = zeros(I,code_length);
    for i = 1:I
        for m = 1:code_length
            R(i,m) = G(i,m) * code (n,m);
            tot(i) = tot(i) + R (i,m);
         end
    end
    data_received = [data_received ; tot / code_length];
end
disp('Reconstructed data at the receiver = '),disp(data_received);

Output:

Enter Code Length = 4
Enter Data Bit Stream = [ 1 -1  1 -1  1  1 -1 -1 ;
        -1 -1  1  1  1 -1 -1  1 ;
         1  1 -1 -1 -1  1  1 -1 ;
         1  1  1  1 -1 -1 -1 -1 ];

Walsh code generated =
    -1    -1    -1    -1
    -1     1    -1     1
    -1    -1     1     1
    -1     1     1    -1

Data bit stream to be transmitted =
     1    -1     1    -1     1     1    -1    -1
    -1    -1     1     1     1    -1    -1     1
     1     1    -1    -1    -1     1     1    -1
     1     1     1     1    -1    -1    -1    -1

Walsh codes used for transmission =
    -1    -1    -1    -1
    -1     1    -1     1
    -1    -1     1     1
    -1     1     1    -1

MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm


texts: Why is peak detection in chromatography not completely automatic?

texts:
When I started to write a software for chromatogram processing I naively thought I'd find some existing algorithms that do all the peak detection automatically and there won't be a need for any manual curation whatsoever. But it appeared that none of the major software in this field works this way. Users have to manually set up parameters (peak widths, thresholds) which work on some chromatograms but not others. So it's either an endless tuning of parameters or manual correction.
My question is philosophical: what keeps this problem unsolved? Do other fields (astronomy, particle physics) have similar problems interpreting measurements automatically? Can other approaches except DSP (e.g. Machine Learning) be useful here?


texts:
When creating a scipy chirp function I get some clicks at the start end the end.
The first one is quite easy to fix using a phase difference of 90 degrees.

The problem I have is that the end point is not 0.

for an static signal I can use the following code:

period_samples = fs / f
periods = np.ceil(T * fs / period_samples)
T = periods * period_samples / fs

and then I have the next zero point where I can stop.

But how can I do something similar with a chirp

my code now is:

import numpy as np
from scipy.signal import chirp

T = 1  # sec
fs = 44100  # samplerate
f0 = 20  # Hz
f1 = 20000  # Hz

t = np.arange(T * fs) / fs
Sig = chirp(t, f0, T, f1, 'logarithmic', 90)


texts:
I have two small radar sensors that use 25 MHz. From each sensor I get the sent signal and received signal (Doppler effect). I know that I can use the FFT, use a window, or FIR filter, take the inverse, and get the frequency within that window. Thus I can analyze the signals separately. I'll have two plots. The radars can be positioned in 3D space at different points, but they can be pointed at the same point in space or at different points.

Is there any way to combine the two signals to analyze motion? For example if one sensor is pointed at the cheek and the second sensor is pointed at the chin, can I somehow average the signals to get one signal that represents the observation of the cheek and the chin? Can I average the two separate plots together?

Or say we had a spinning top. One sensor is on the right and another on the left. If the top leans right the sensor on the right gets less of a time delay and the sensor on the left gets a longer time delay for the received signal. Thus using the information from both I can tell at what angle the top is spinning hypothetically because speed is distance over time.

Any help, math equations or references would be appreciated.


texts: Deriving expression for the DTFT of a rectangular window

texts:
Looking at the picture above, how did the author get from point A) to B)?

My Approach: Multiply A) by $e^{j\omega/2}/e^{j\omega/2}$. Now I am stuck with simplying the numerator.


texts:
I have implemented FIR filter using tapped delay line method. I start getting output as soon as first input sample is passed to it, I am wondering from which sample I will get the proper output(without transient outputs) from filter.

Does it depend on location of highest magnitude tap or it's just dependent on the length of filter.The number of filter coefficients I am using are N(odd) and symmetric.
for verifying I took the following filter taps as example :-

$F1$ $=$ $[1 0 0 0 0 0 0];$

$F2$ $=$ $[0 1 0 0 0 0 0];$

$F3$ $=$ $[0 0 1 0 0 0 0];$
and so on,if I pass sine wave using filter coefficients $F1$ then the first sample of the output is itself meaningful but if I use $F2$ and $F3$ then meaningful output is from second and third output samples respectively. I am also wondering how matlab command $filter()$ works.


texts:
I have DSP in my academics , I am here to have a suggestion regarding the better alternative to MATLAB because of Its downloading size.

Thank you.


texts:
http://zone.ni.com/images/reference/en-XX/help/371025N-01/loc_eps_rectangular_modulation.gif

I just wanted to make sure that my assumption that in this figure $\phi(t)$ is representing $m(t)$, and $R(t)$ is  the amplitude of the carrier frequency is correct.


texts: Why sinc function defined on minus side?

texts:
What this relation on picture mean? Why response for a digital low-pass filter exists on "-" side of x axis? Impulse starts at 0 as i can understand and goes to infinity...
Please use phrases starting with "think of it this way..." or "it's like..."?


texts:
I am somewhat confused on the terminology used in Computer Vision. Suppose I wanted to detect (localize) an object (say, a human) in an image. However, a box around the human is not sufficient, what I need is a pixelwise localization. More precisely, I need a mask telling me which pixels belong to the human and which ones don't.
If I wanted to find papers dealing with this problem, what would be the term to search for?

Object detection: Yes, but most approaches actually work in terms of bounding boxes
Object localization: If I see it correctly, this is essentially used as a synonym for object detection
Image segmentation: Here we segment objects in terms of pixels but we do not necessarily only aim to segment out one specific object that is known.

Is there very famous work in this area I should not miss?


texts: Which IIR filters approximate a Gaussian filter?

texts:
So it recently dawned on me that Bessel filters, despite being listed along with the other common types, are really an oddball that belongs in a different "class", and I'm trying to learn more about it.

  The rectangular magnitude response represents the ideal frequency domain response, for the transition band is zero and the stopband has infinite attenuation. The Gaussian magnitude response, on the other hand, represents the ideal time-domain response, in that no overshoots occur in the impulse response and the step response. Many of the responses attained in practice are approximations to these ideal ones  source

So a brickwall filter is convolution with a sinc function, and has these frequency domain properties:

Flat passband
Zero stopband
Infinite roll-off rate/no transition band

It's non-causal and unrealizable because of the infinite tails in both directions.  It is approximated by these IIR filters, with the approximation improving as order increases:

Butterworth (maximally flat passband)
Chebyshev (maximum roll-off rate with stopband or passband ripple)
Elliptic (maximum roll-off rate with stopband and passband ripple)
Legendre  (maximum roll-off rate with monotonic passband)

The Gaussian filter is convolution with a Gaussian function, and has these time domain properties:

Zero overshoot
Minimal rise and fall time
Minimal group delay

It's unrealizable for the same reasons as the sinc function, and can be approximated by these IIR filters, more closely as order increases:

Bessel (maximally flat group delay)  according to 1 and 2

Here's Bessel filters of increasing order along with a Gaussian dashed line I picked merely because it seemed to fit the trend ($e^{-{1 \over 2}(\pi \omega)^2}$):

So my questions are:

Is everything right so far?  If so, are there other IIR filters that approximate the Gaussian?  What are they optimized for?  Maybe one that minimizes overshoot?

If you search for "IIR Gaussian" you can find a few things (Deriche? van Vliet?), but I don't know if they're really the same as a Bessel or if they optimize for some other property, etc.


texts:
Setup

Suppose we have a complex $L\times 1$ signal $\mathbf{x}$ with two tones at (unknown) frequencies and phases defined as:
$$
x_n = A_1 e^{j \omega_1n + \varphi_1} + A_2 e^{j \omega_2n + \varphi_2}
$$
for $0 \leq n \leq L-1$. We observe $\mathbf{x}$ corrupted by additive white Gaussian noise:
$$
\mathbf{y} = \mathbf{x} + \mathbf{w}
$$
where $\mathbf{w} \sim \mathcal{CN}(0, \sigma^2 \mathbf{I})$.

Problem

I'm trying to show a connection between the signal amplitudes and eigenvalues from the eigendecomposition
of the autocovariance matrix $\mathbf{R_y}$ of $\mathbf{y}$ (this is related to MUltiple SIgnal Classification i.e. MUSIC algorithm which splits the $L$ dimensional space into a 2 dimensional signal subspace and an $L-2$ dimensional noise subspace).

$$
\mathbf{R_y} := \mathbf{E}[\mathbf{y}\mathbf{y}^H] = \mathbf{x}\mathbf{x}^H + \sigma^2 \mathbf{I}.
$$

Let $\mathbf{s}_1 = \begin{bmatrix} 1 & e^{j\omega_1} & e^{j2\omega_1} & \ldots & e^{j(L-1)\omega_1}\end{bmatrix}^T$
and similarly define $\mathbf{s}_2$. Then,
$$
\mathbf{x} = \begin{bmatrix}
           \mathbf{s}_1 & \mathbf{s}_2
         \end{bmatrix}
         \begin{bmatrix}
         A_1 e^{j\varphi_1} \\
         A_2 e^{j\varphi_2}
         \end{bmatrix}
$$
and the signal autocovariance matrix can be written as:
\begin{eqnarray*}
\mathbf{x}\mathbf{x}^H &=& \begin{bmatrix}
           \mathbf{s}_1 & \mathbf{s}_2
         \end{bmatrix}
         \begin{bmatrix}
         A_1 e^{j\varphi_1} \\
         A_2 e^{j\varphi_2}
         \end{bmatrix}
         \begin{bmatrix}
         A_1 e^{-j\varphi_1} & A_2 e^{-j\varphi_2}
         \end{bmatrix}
         \begin{bmatrix}
           \mathbf{s}_1^H \\
           \mathbf{s}_2^H
         \end{bmatrix} \\
         &=& \underbrace{\begin{bmatrix}
           \mathbf{s}_1 & \mathbf{s}_2
         \end{bmatrix}}_{\mathbf{S}}
         \underbrace{\begin{bmatrix}
         A_1^2 & A_1 A_2 e^{j(\varphi_1-\varphi_2)} \\
         A_1 A_2 e^{-j(\varphi_1-\varphi_2)} & A_2^2
  \end{bmatrix}}_{\mathbf{M}}
  \underbrace{\begin{bmatrix}
           \mathbf{s}_1^H \\
           \mathbf{s}_2^H
         \end{bmatrix}}_{\mathbf{S}^H}
\end{eqnarray*}
where we note that $\mathbf{M}$ is Hermitian symmetric and therefore has a decomposition $\mathbf{M} = \mathbf{Q\Lambda Q}^H$ where $\mathbf{\Lambda}$ is
a 2x2 diagonal matrix and $\mathbf{Q}^H\mathbf{Q} = \mathbf{Q}\mathbf{Q}^H = \mathbf{I}$.

Therefore,
\begin{eqnarray*}
\mathbf{R_y} &=& \mathbf{SQ \Lambda Q^HS^H} + \sigma^2 \mathbf{I}.
\end{eqnarray*}

Now I am stuck. I suspect I should be able to eventually express $\mathbf{R_y} = \mathbf{PDP^H}$ where the diagonal matrix $\mathbf{D}$ should have two large eigenvalues related to the signal amplitudes $A_1$ and $A_2$, but since $\mathbf{SS^H}\neq \mathbf{I}$ I have no way to absorb the $\sigma^2$'s into the diagonal matrix $\mathbf{\Lambda}$.


texts: Maximally Stable Extremal Regions (MSER) Implementation in Document Image Character Patch Identification

texts:
My task is to identify character patches within the document image. Consider the image below:

Based from the paper titled "CG-DIQA: No-reference Document Image Quality Assessment Based on Character Gradient", to extract character patches, the MSER based method will be adopted to detect character candidates.

  "The main advantage of the MSER based method is that such algorithm is
  able to find most legible characters even when the document image is in
  low quality."

Another paper titled "Robust wide-baseline stereo from maximally stable extremal regions" discusses about MSER. I'm having a hard time understanding the latter paper. Can anyone explain to me in simple terms the steps that I should take to implement MSER and extract character patches in my sample document. I will implement it in Python, and I need to fully grasp / understand how MSER works.

Below are the steps to identify character patches in the image document (based from the way I understand it, please correct me if I am wrong)

"First, pixels are sorted by intensity"

My comprehension:

Say for example I have 5 pixels in an image with intensities (Pixel 1)1, (Pixel 2)9,(Pixel 3)255,(Pixel 4)3,(Pixel 5)4 consecutively, then if sorted increasingly, based on intensity it will yield an output, Pixel 1,4,5,2 and 3.
After sorting, pixels are placed in the image (either in decreasing or increasing order) and the list of connected components and their areas is maintained using the efficient union-find algorithm.

My Comprehension:

Using the example in number 1. Pixels will be arranged like below. Pixel component/group and Image X,Y coordinates are just examples.

 Pixel Number | Intensity Level | Pixel Component/Group | Image X,Y Coordinates
      1       |        1        |  Pixel Component # 5 | (14,12)
      4       |        3        |  Pixel Component # 1 | (234,213)
      5       |        4        |  Pixel Component # 2 | (231,14)
      2       |        9        |  Pixel Component # 3 | (23,21)
      3       |      255        |  Pixel Component # 1 | (234,214)

"The process produces a data structure storing the area of each connected component as a function of intensity."

My comprehension:

A column in table in #2 will be added called Area. It will count the number of pixels in a specific component with the same intensity level.

4."Finally, intensity levels that are local minima of the rate of change of the area function are selected as thresholds producing MSER. In the output, each MSER is represented by position of a local intensity minimum (or maximum) and a threshold."

Questions:

How to get the local minima  of the rate of change of the area function ?

Please help me understand this.


texts: Kalman Filter : How measurement noise covariance matrix and process noise helps in working of kalman filter , explain intuitively please?

texts: How noise covariance matrix and process covariance matrix helps in improving the state estimate, can some one explain intuitively without mathematics ?


texts: The Rank Of Broadband LCMV Beamformer Constraint Matrix

texts:
So I am wondering why I keep seeing that the rank of the LCMV Broadband Constraint Matrix is equal to the number of constraints, when in practice, I find that the rank is only 1. This is precluding me from performing any sort of 'low rank' approximation since it is already of the lowest rank possible!

The constraint matrix is defined as:

$C = [ \mathbf{d}(\theta, \omega_{0}) | \mathbf{d}(\theta, \omega_{1})| \ldots |\mathbf{d}(\theta, \omega_{r-1})] $

Where each $\mathbf{d}(\theta, \omega_{i})$ is a column vector of the form:

$ \mathbf{d}(\theta, \omega_{i}) = [ e^{-j \omega_{i} \tau_{0}}, \ldots, e^{-j \omega_{i} \tau_{M-1}} | e^{-j \omega_{i} (\tau_{0} + T_{s})}, \ldots, e^{-j \omega_{i} (\tau_{M-1} + T_{s})} | \ldots | e^{-j \omega_{i} (\tau_{0} + (J-1)T_{s})}, \ldots, e^{-j \omega_{i} (\tau_{M-1} + (J-1)T_{s})}]^{T}$

Where $M$ is the number of sensors, $J$ is the length of the attached FIR filter, $T_{s}$ is the sampling interval, and $\tau_{i}$ is the delay of the signal coming from angle $\theta$ to the $i$th sensor.

When I create this matrix and compute the SVD via zgesvd_ , I get only one significant singular value. The literature says the rank should be $r$, ie the same as the number of columns of $C$, but this is not the case -- ie the columns have linear dependence.

Why do they say this? Is there something I am missing? Has anyone else encountered this?

References:

Wei Liu and Stephan Weiss. Wideband Beamforming: Concepts and Techniques. Wiley. 2010. pp. 33-35

Kevin M. Buckley. Spatial/Spectral Filtering with Linearly Constrained Minimum Variance Beamformers. IEEE Trans. Acoustics, Speech, and Signal Processing. March 1987 pp.249-266


texts: Finding polynomial approximations of a sine wave

texts:
I want to  approximate the sine wave given by $\sin\left(\pi x\right)$ by applying a polynomial waveshaper to a simple triangle wave, generated by the function

$$T\left(x\right)=1-4\left|\tfrac{1}{2}-\operatorname{mod}(\tfrac{1}{2}x+\tfrac{1}{4},\ 1)\right|$$

where $\operatorname{mod}(x, 1)$ is the fractional part of $x$:

$$ \operatorname{mod}(x, y) \triangleq  y \cdot \left( \left\lfloor \frac{x}{y}\right\rfloor - \frac{x}{y} \right) $$

A Taylor series could be used as a waveshaper.

$$S_1\left(x\right)=\frac{\pi x}{2}-\frac{\frac{\pi x}{2}^3}{3!}+\frac{\frac{\pi x}{2}^5}{5!}-\frac{\frac{\pi x}{2}^7}{7!}$$

Given the above functions, $S_1(T(x))$ will get us a decent approximation of a sine wave. But we need to go up to the 7th power of the series to get a reasonable result, and the peaks are a little low and won't have a slope of exactly zero.

Instead of the Taylor series, we could use a polynomial waveshaper following a few rules.

Must pass through -1,-1 and +1,+1.
Slope at -1,-1 and +1,+1 must be zero.
Must be symmetrical.

A simple function meeting our requirements:

$$S_2\left(x\right)=\frac{3x}{2}-\frac{x^3}{2}$$

The graphs of $S_2(T(x))$ and $\sin\left(\pi x\right)$ are pretty close, but not as close as the Taylor series. In between the peaks and zero-crossings they visibly deviate a bit. A heavier and more accurate function meeting our requirements:

$$S_3\left(x\right)=\frac{x(x^2-5)^2}{16}$$

This is probably close enough for my purposes, but I'm left wondering if another function exists which approximates the sine wave more closely, and is computationally cheaper. I have a pretty good grasp on how to find functions meeting the three requirements above, but I'm not sure how go about finding functions which meet those requirements and also most closely match a sine wave.

What methods exist for finding polynomials which mimic a sine wave (when applied to a triangle wave)?

To clarify, I'm not necessarily looking only for odd-symmetrical polynomials, although those are the most straightforward choice.

Something like the following function could also fit my needs:

$$S_4\left(x\right)=\frac{3x}{2}+\frac{x^2}{4}+\frac{x^4}{4}$$

This meets the requirements in the negative range, and a piecewise solution could be used to apply it to the positive range as well; for example

$$\frac{3x}{2}-\frac{P\left(x,2\right)}{4}-\frac{P\left(x,4\right)}{4}$$

where $P$ is the signed power function.

I'd also be interested in solutions using the signed power function to support fractional exponents, as this gives us another "knob to twist" without adding another coefficient.

$$a_0x\ +a_1P\left(x,\ p_1\right)$$

Given the right constants, this could potentially get very good accuracy without the heaviness of fifth- or seventh-order polynomials. Here's an example meeting the requirements described here using some hand-picked constants: $a_0=1.\overline{666}, a_1=-0.\overline{666}, p_1=2.5$.

$$\frac{5x-2P\left(x,\ \frac{5}{2}\right)}{3}$$

In fact, those constants are very close to $\frac \pi 2$, and $1 - \frac \pi 2$, and $e$. Plugging those in gives something that looks extremely close to a sine wave.

$$\frac{\pi}{2}x\ +\left(1-\frac{\pi}{2}\right)P\left(x,e\right)$$

To put it another way, $x-\frac{x^e}{6}$ looks very close to $\sin(x)$ between 0,0 and π/2,1.  Any thoughts on the significance of this? Maybe a tool like Octave can help discover the "best" constants for this approach.


texts: Sampling rate versus number of sample of I and Q

texts:
I am acquiring data of an O-QPSK signal (chip rate : 2Mchip/s) with a SDR. I am undersampling the signal at 20 Mhz (carrier is at 2.8 GHz) so I respect the Shannon-Nyquist criterion that you need to sample 2 time the bandwidth of your signal, however when I take apart the I and Q phase, the signal look like it has been sampled at 10 Mhz. I observe this when I create my reference signal in matlab with a sampling rate of 20 Mhz for each phase I&Q. I end up with a signal x2 time the size of the real signal (2560 samples vs 5120 samples).

The chip period of my signal is 1 microsecond and I am sampling at a rate of a rate of 20 Mhz  with 128 in my sampling period, you effectively end up with 2560 samples.

When you take each phase alone, each part of a chip will be 2 time the period chip which result in a signal of 2 microsecond, therefore is we don't divide the sampling rate by 2, we will end with 5120 sample.

This effect seem to be produce by the phase offset that the O-QPSK use. Can someone give an explanation about this phenomenon?

[EDIT]

The SDR I am using is the B200 Mini from Ettus Research.


texts:
I'm doing a project where I need to measure the rate of vibrato cycles in a sung note under various conditions. The stimulus will be a several-second-long audio sample of a held a capella note with vibrato; the output should be the rate of the vibrato repetition (e.g. 4Hz, 10.5Hz, etc). The audio will be captured specifically for this study; I won't be sampling existing audio, and don't need to deal with interfering signals. The underlying pitch of the note can be whatever is needed, although precision probably won't be possible.

Here's an example sonogram:

Ideally, this would be done with a Linux command line chain or shell script, but if I need to write some Python that would be fine, too.


texts: Histogram of signal with a time shift added

texts:
One of the problems in Proakis's signal processing book using Matlab asks to plot the histogram of

x3(n) = x1(n) + x1(n-1)

where x1(n) is a random sequence.

To create x1(n), I do the following:

n1 = 1:100000;
a = 0; b = 2; % interval over which to pick random numbers.
x1 = a + (b-a)*rand(100000,1);
N1 = hist(x1,100);
figure, bar(N1);

To generate x3(n), I use a couple of helper-functions defined in Proakis's book.

[x31,n31] = sigshift(x1,n1,1);
[x3,n3] = sigadd(x1,n1,x31,n31);
N3 = hist(x3,nbins);
figure, bar(N3);

Unfortunately, I do not have enough points to post images. However, the histogram of the random sequence x1(n) just shows that x1 is uniformly distributed. The x3(n) however, looks more like a Gaussian distribution.

I don't understand why this happens, and would appreciate if someone could explain this a bit.

Thanks.


texts: Color Artifacts in Fourier Transformed Image

texts:
I am working with images and Fourier transforms.  I am trying to understand what might be causing some artifacts in my output image.

I am starting with a 512x512, RGB image of Lenna.  I FFT the image, and run a low-pass filter on it.  Unchanged for everything inside the circle, 0 for everything outside the circle.  Here is what the mask looks like:

So far so good.  When I inverse FFT back to an image, I get bright green artifacts:

These look out of place to me, what might be causing them?  Is it the sharp cutoff in my filter?  I would have guessed that would creating ringing artifacts, but nothing like a sharp color blob.  Also why green only?  It's not like there's a lot of green in the original.

I'm using numpy on python 2.7.6 in 64-bit Linux, if that matters.


texts: Picking the correct filter for accelerometer data

texts:
I am fairly new to DSP, and have done some research on possible filters for smoothing accelerometer data in python. An example of the type of data Ill be experiencing can be seen in the following image:

Essentially, I am looking for advice as to smooth this data to eventually convert it into velocity and displacement. I understand that accelerometers from mobile phones are extremely noisy.

I dont think I can use a Kalman filter at the moment because I cant get hold of the device to reference the noise produced by the data (I read that its essential to place the device flat and find the amount of noise from those readings?)

FFT has produced some interesting results. One of my attempts was to FFT the acceleration signal, then render low frequencies to have a absolute FFT value of 0. Then I used omega arithmetic and inverse FFT to gain a plot for velocity. The results were as follows:

Is this a good way to go about things? I am trying to remove the overall noisy nature of the signal but obvious peaks such as at around 80 seconds need to be identified.

I have also tired using a low pass filter on the original accelerometer data, which has done a great job of smoothing it, but I'm not really sure where to go from here.
Any guidance on where to go from here would be really helpful!

EDIT: A little bit of code:

for i in range(len(fz)):
    testing = (abs(Sz[i]))/Nz

    if fz[i] < 0.05:
        Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)
Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

So essentially, ive performed a FFT on my accelerometer data, giving Sz, filtered high frequencies out using a simple brick wall filter (I know its not ideal). Then ive use omega arithmetic on the FFT of the data.
Also thanks very much to datageist for adding my images into my post :)


texts: Adaptive volume normalization on a stream of PCM data

texts:
I am writing a program that receives a continual stream of PCM data from various sources. Some sources have different volumes than others. The PCM data is always a live continuous stream of 16 not signed data single channel.

How can I dynamically make sure that all samples are not too low in volume or high in volume to effectively normalize its output realtime?


texts: Estimation of step amplitude using the CUSUM algorithm

texts:
I used the CUSUM algorithm to detect steps in data. Basically the data looks like this, the data has a constant amplitude and then there is a rapid variation or a step. For example, the signal has a value of 1000 and rapidly drops to 950 or 900. With the CUSUM algorithm, I can detect easily these changes and the time at which they occur. But now I would also like to automate the amplitude estimations of these steps. Eventually, I'll need to implement this in real-time ( a small delay is tolerable for the amplitude estimation)

One solution I would like to try is this one :

take N samples before the step and take N samples after the step and subtract the respective averages.

Is there a better solution?

Edit :
- The transition from 1000 to 900, for example, might take 10 samples in all.


texts:
I am trying to understand this equation:

It comes from:

$Force = mass * acceleration$

$F(t) = m * a(t)$

$F(s) = m * (s^2 * y(s) - s*y0 - y0)$, where $y0=0$

$F(s) = m * s^2 * y(s)$

$y(s) =  F(s)/(m*s^2)$

Then I don't understand what type of substitution they are using for s. It appears they substitute:

$s =  (1-z^{-1})/T$

And then the last $z^{-1}$ is just to indicate they are delaying the whole thing one sample for another reason.

What type of s substitution is this called?


texts: Perform inverse distance weighting interpolation using multiple images matlab

texts:
I have a total of 4 highly identical grayscale images, where image1 is taken as the reference image, the rest of the images are to be used to enhance image1.

My objective is to enlarge the image to 2 times the original size. Instead of taking average values from neighbouring pixels to fill the upsampled pixels, is it possible to do inverse distance weighting interpolation, so that my image remains sharp?

My other concern is that I need to fill the upsampled pixels with the pixel values from the 3 other images. Is it possible to do interpolation with 4 input images?


texts:
I'm really tripping about this point! I try to reason with you and explain what's my trouble.
Consider a scenario with Sample Rate 44100hz and a LFO square wave that will modulate some internal parameter of my plugin.

Its clear that if I generate that signal using this function:

if (mPhase <= mPI) {
    value = 1.0;
}
else {
    value = -1.0;
}

mPhase += mPhaseIncrement;
while (mPhase >= mTwoPI) {
   mPhase -= mTwoPI;
}

once it reaches my D/A converter, the interpolation would make a signal that will contains alias (the more I pitch up foundamental, the more alias I got).

That's normal: I don't have sufficient samples to represent all partials of the signal I'd like to represent. Clear.

But, why this is true also on LFO?

I mean: once I create my digital square (that, until it interpolate, it's really just 1/-1), the target param would be modulated really as a "real" square wave.

Why it should modulate the knob in a way that contains "alias"? The knob won't move as the signal being generated after the interpolation (with ripples and not linear movement). It really moves linear as a square wave: -1, 1, -1, 1. No "alias".

I hope you got what I mean? I can't see visually the "alias" introduced by LFO (by square wave, in this case) modulating a knob.

Note: I know that using an LFO with the code above at a frequency lower than 20hz produce "un-noticeable" alias for audio/human hearing. But for the reason I question, it shouldn't produce alias at all.


texts: Notation confusion -- What is the correct operator for computation of the log-likelihood expression for complex valued data?

texts:
This question is an extension of another question of mine asked earlier here Help in understanding if the maximum likelhood estimation is working properly

In that question the inputs were real valued. Now, in this question the inputs are real valued and hence the confusion.

I am having difficulties in the notation when applying to the complex domain.
What notations and operators should I use when all the data are complex valued? Please help


texts: Shannon capacity with distribution different from Gaussian

texts:
If I understand correctly, in AWGN channel, for a given SNR, there is a code for $M$ input codeword $\mathbf{X}_{i, 1 \leq i \leq M}$ that
$$\lim_{n \to \infty} R=\lim_{n \to \infty}\frac{\log(M)}{n} =\log(1+SNR)$$

(where $n$ is the length of $\mathbf{X}_i$) if $\mathbf{X}_i$ is drawn from multivariate Normal distribution.

However, I don't know if $\mathbf{X}_i$ is not Normal, does the Shannon channel coding theorem hold ?


texts:
Parallel programming is often used to accelerate image processing methods (e.g. Eklund et al., 2013); however, there doesn't appear to be the same level of discussion about using parallel processing to handle gigapixel images (or images that are too large to load on one machine). Is there a reference that discusses parallel processing techniques for large images?

Follow-up Question: To me, it seems there is only one way to process gigapixel images: break them up into chunks and process each chunk on a different node. The way you split the image could differ, but the basic procedure stays the same. Are there other ways parallel processing can be used to handle gigapixel images?


texts:
In below image, we have scaling property of DFT, how the final equation is obtained from the above equation. That is how we are getting the scaling factor , $ \frac{1}{|ab|} $ in the final equation ? Please explain how we are getting $ \frac{1}{|ab|} $. Thank you

Equation 4.12 is given as


texts: RMS typical length for audio dynamic compressor

texts:
I know what does it mean Root Mean Square. But I can’t find any information what is typical portion (length) of audio signal for RMS envelope detector in dynamic audio compressor.

Of course I know it depends on compressor design, there are many compressors with RMS detection. But I wonder is there any reccomendation what gives resonable good results? Like for example should it be longer than 1 second or less?
I mean some compression for audio mastering or mix busses?

I don’t have any idea, should it be just little bit of second like 0.1 second? Or better more seconds like 2.0 or even 3.0 seconds?

Could anyone give some advice? I would be very appretiate.

And what is also interesting, how about attack and release time for RMS detector compression? It seems to be stupid to use shorter attack time than length of RMS portion pf audio. But maybe it is OK?


texts: Savitzky-Golay smoothing filter for not equally spaced data

texts:
I have a signal that is measured at 100Hz and I need to apply the Savitzky-Golay smoothing filter on this signal. However, on closer inspection my signal is not measured at perfectly constant rate, the delta between measurements ranges between 9.7 and 10.3 ms.

Is there a way to use the Savitzky-Golay filter on not equally spaced data?
Are there other methods that I could apply?


texts: Computing the decimation ratio between two m-sequences

texts:
Let's suppose I have an LFSR that generates an m-sequence $y_1[k]$ --- in other words, the LFSR has $N$ bits and $y_1[k]$ has period $m=2^N - 1$.

Now suppose I know someone has decimated this and taken every $j$th element, so $y_2[k] = y_1[jk+b]$. And I know the $y_2[k]$ but the values of $j$ and $b$ are unknown.

Is there any efficient way to compute the decimation ratio $j$ without testing all possible values of $j$?

(Note 1: if $b$ is unknown then there are $N$ possible solutions for $j$ because $2j, 4j, 8j, \ldots\, j^{2^{N-1}}$ can have sequences generated by LFSRs with the same coefficients, differing only by a time shift. I have no idea if the problem is easier to solve for known $b$.)

(Note 2: This smells vaguely like a discrete logarithm problem but I cannot figure out how to restate it as such. I will be satisfied if I can obtain an equivalent discrete logarithm problem, because the values of $N$ that I deal with are small enough that I can use established techniques like Pollard's rho and kangaroo algorithms to compute the value of $j$.)

There are several equivalent ways of stating this problem (I posted one on https://math.stackexchange.com/questions/2480132/determining-decimation-ratio-given-characteristic-polynomials-of-quotient-rings):

$y_1[k] = \operatorname{Tr}_{F_1}(w_1x^k)$ where $F_1 = GF(2)[x]/p_1(x)$ and similarly for $y_2[k], F_2, w_2, p_2(x)$; the traces are related (because of the decimation ratio $j$) by $\operatorname{Tr}_{F_1}(w_1x^{jk+b}) = \operatorname{Tr}_{F_2}(w_2x^k)$; the coefficients of $p_1(x)$ are the LFSR coefficients; $p_2(x)$ can be determined from the Berlekamp-Massey algorithm, and $w_1$ and $w_2$ can be determined from state recovery methods by running the LFSR backwards (I have a blog article in progress on this but it's unfinished and not a quick thing to summarize, sorry)
$p_1(x)$ is the minimal polynomial of $x$ in $F_1 = GF(2)[x]/p_1(x)$; $p_2(x)$ is the minimal polynomial of $x^j$ in $F_1$ and $p_2(x)$ is also the minimal polynomial of $x$ in $F_2 = GF(2)[x]/p_2(x)$; given $p_1(x)$ and $p_2(x)$, what is $j$?

Example:

$p_1(x) = x^9 + x^4 + 1$
$y_1[0..18] = [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1]$
$j = 107, b = 80$
$y_2[0..18] = [1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0]$

We can use Berlekamp-Massey on the $y_2$ values to determine $p_2(x) = x^9+x^7+x^5+x+1$, and state recovery techniques to determine that $y_2[k] = HB(v_2x^k)$ with $v_2 = x^8 + x^6 + x^5 + x^4 + x^3$ (equivalently the LFSR has initial state 101111000), where $HB(u)$ is the high bit = coefficient of $x^8$.

But I'm stuck on how to get from there to $j=107$ or other members of its cyclotomic coset $\{107, 214, 428, 345, 179, 358, 205, 410, 309\}$, aside from checking every value of $j$ --- or at least, the smallest coset representative of every value of $j$, and then using Berlekamp-Massey until I get a match of $p_2(x)$.


texts: Why do AAC files seem to have an inexact bitrate?

texts:
I'm using ffmpeg to encode audio files as AAC at 256kbps. However, the actual bitrate of the resulting file is rarely 256. It's usually a number close to 256 like 241, 255, 257, 254. In one case it was 125.

What's going on? Is this is a property of AAC or ffmpeg? If I ask ffmpeg for mp3 at 320kbps I always get 320. I thought for a while that the AAC encoder was using 256 as an upper limit and only using the bitrate it needs but then I got a few 257 files so that can't be it.


texts: What is the effect of quantization bucket width in Locality Sensitive Hashing?

texts:
In the LSH Tutorial by Slaney, it is written that

  Increasing the quantization bucket width w will increase the number of
  points that fall into each bucket. To obtain our final nearest
  neighbor result we will have to perform a linear search through all
  the points that fall into the same bucket as the query, so varying w
  effects a trade-off between a larger table with a smaller final linear
  search, or a more compact table with more points to consider in the
  final search.

As per my understanding, if the bucket size is large, then there should be less number of tables (each bucket is a table right ? ). Since, each point is hashed to a cuket. If a query is hashed to a bucket, then with a linear scan on a larger table, wouldn't it encounter more points than a compact table ?


texts: How to treat noise (in the acceleration) in a Kalman filter when tracking position?

texts:
I have an equation of motion for which I am wanting to track the position.

The equation looks like this:

$$\ddot{x}(t) = -a\dot{x}(t) - bx(t) + F(t)$$

where $F(t)$ is an external Stochastic force noise (like brownian motion).

I know if $F = 0$ I can build a Kalman filter like so:

$$
\begin{bmatrix}
\dot{x} \\
\ddot{x} \\
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
-b & -a \\
\end{bmatrix}
.
\begin{bmatrix}
x \\
\dot{x} \\
\end{bmatrix}
$$

And then transform $ \begin{bmatrix}
0 & 1 \\
-b & -a \\
\end{bmatrix} $ like so:

\begin{align}
A(t) & = \mathscr{L}^{-1}\left\{\left(s\vec{I} - \begin{bmatrix}
0 & 1 \\
-b & -a \\ \end{bmatrix}\right)^{-1}\right\}\\
& = \mathscr{L}^{-1}\left\{\left(\begin{bmatrix}
s & 0 \\
0 & s \\ \end{bmatrix} - \begin{bmatrix}
0 & 1 \\
-b & -a \\ \end{bmatrix}\right)^{-1}\right\}\\
&= \mathscr{L}^{-1}\left\{\left(\begin{bmatrix}
s & -1 \\
b & s+a \\ \end{bmatrix}\right)^{-1}\right\}
\end{align}

... (at this point I get stuck as I am using Sympy in python to calculate the inverse Laplace transform and it fails for this matrix).

But I want to know how I should treat this Stochastic effect in a Kalman filter.


texts: variable length no-refrence filter algorithm

texts:
I'm searching for a an algorithm or an idea of implementation. I need to filter a signal in real time, and want to change the length of it (number of taps) in order to change the delay.

In low signal environment I want a high order filter and care less
about delay
In high signal environment I want a low order filter in order to
minimize the delay.

I've tried to search "variable length filter..." on google but all of the results were versions of the LMS algorithm, but I don't have a "clean noise" reference signal which I can use.

Have anyone heard of something like this?
Thanks :)


texts: How do I find the Hilbert transform (qudrature) of a discrete signal by convolution with a Hilbert kernel?

texts:
Given a discrete real signal $ f_{re}(t) $ the (complex) analytic signal is given by:
$ x(t) = f_{re}(t) + j*f_{im}(t) $.
I want to calculate $f_{im}(t)$: the quadrature by convolving with a Hilbert kernel of size 2*n+1:
$f_{im}(t) = H*f_{re}(t)$
$H=[H_{-n},H_{-n+1},..,H_0,H_{n-1},H_n]$
According to "Complex seismic trace analysis",M. T. Taner*, F. Koehler*, and R. E. Sheriff;
$ H_i = 2*sin^2(\pi i/2)/(\pi i), H_0 = 0$
for $ n = \infty $.
I want to use a finite length kernel, e.g. n=16.
Do I need to taper the kernel from the equation above with a window?
How do I normalize the odd (antisymmetric) result kernel?
Am I right in assuming that calculating the quadrature by convolving with H introduces a smoothing?
If so; should I also smooth the real signal, e.g. by a gaussian kernel of size 2n+1?
EDIT:
In Matlab there is a function called "firpm" that creates a "Parks-McClellan optimal equiripple FIR filter".
It seems this is generally considered the "optimal" FIR approximation to an IIR?
Does this mean that using the simple coefficient equation from "Complex seismic trace analysis" is too simplistic?
EDIT2:
regarding the normalization: I tested in Matlab and found that convolving with a kernel calculated from the simple coefficient equation gave good results in the case of the input signal being a sine. However I had to divide by a factor m to get correct normalization:
H5 = [-0.0000   -0.6366  0    0.6366    0.0000]: m = 0.1271
H7 = [-0.2122 -0.0000 -0.6366  0  0.6366 0.0000 0.2122]: m=0.2525

What is the general formulae for m?
EDIT 3:
The ideal kernel would have constant magnitude response. However it seems that since my kernel will be of finite size; it will act as an bandpass filter, and therefore I need to look out for bandpass ripple.
The solution seems to be a magnitude response that rolls of gently at the low and high frequencies. The "simple coefficient equation" is quite good in this respect; it does not introduce a lot of ripple, instead the shorter the kernel the earlier it starts rolling off. As an example here is the 19 point kernel:

In comparison the 19 point kernel from Matlabs firpm does not roll of at high frequencies (!?) but otherwise have roughly the same amount of ripple as my simple kernel:

My conclusion so far thus seem to be: no the "simple coefficient equation" gives a good approximation to the optimal FIR for the Hilbert convolution operator. The "Parks-McClellan optimal equiripple FIR filter" might be slightly more "optimal" but the difference is neglible, at least for reasonably large kernels (say 19 points and upwards). Anything I am forgetting? I will therefore move on with this approach. Only remaining problem: how do I normalize?


texts: What is the method for generating coefficients for the Farrow structure in Matlab? I am designing a low pass filter using firpm, but the output of the Farrow structure is not delayed.


texts: What is a periodic signal in image processing?

texts: In the context of image processing (and computer vision), the concept of convolution comes up a lot. Convolution is quite related to the concept of Fourier transform and DFT. In the context of image processing, we are mainly interested in DFTs, given that signals (images) are discrete and finite. In this context, I have heard the expression "periodic signal", which reminds me of a sine or a cosine waves, which have a period of $2\pi$, that is, every $2\pi$ the output of these functions repeats. Now, in the case of images, it is not clear to me the meaning of a periodic signal. What would be a periodic image? How is this concept related to DFT and convolution?


texts:
Do you know of any algorithms for audio processing that make use of peaks in auto-correlation function? I am currently involved in research on mapping methodology for dynamic algorithms. I came recently over an algorithm for detection of fundamental frequency in audio signals, which has a dynamic structure I am looking for - the execution time of the algorithm varies with the number of peaks in the signal's auto-correlation:

while (i < hsize)

    // Search for peaks in auto-correlation:
    if((r[i] > 0) && (r[i]>r[i-1]) && (r[i]>=r[i+1])){

            temp = 0.25f*(r[i+1]-r[i-1]);
            temp = r[i] - temp*temp/(0.5f*(r[i+1]+r[i-1]) - r[i]);
            n = n1[i] + n2[i];
            if ((rmax*n) < (temp*nmax)) {imax=i; rmax=temp; nmax=n;}.

            i++;}
    else{i++;}

This would be perfect as a studying example for my methodology, but the problem is that there is only a couple of lines inside the loop, so the overhead of adding my techniques becomes too big for the problem of this size. Now I am wondering if there exist any similar algorithms in audio processing, i.e. also making use of peaks in auto-correlation function but perform more extensive calculations for each peak?


texts: Is this statement contradicting with the band-limiting theorem?

texts:
The following sentences are from the footnote of page 240 in R.C. Gonzalez's Digital Image Processing (3rd):

  An important special case is when a function that extends from -Inf to
  +Inf is band-limited and periodic. In this case, the function can be truncated and still be band-limited, provided that the truncation
  encompass- es exactly an integral number of periods. A single
  truncated period (and thus the function) can be repre- sented by a set
  of discrete samples satisfying the sampling theorem, taken over the
  truncated interval.

Isn't these contradicting to the band-limiting theorem that a bandlimited signal cannot be also timelimited?


texts:
Newbie question here.

I'm working on a system that is running a 32khz sampling rate, audio, 16 bit DAC.
Previously, our system had a very hard latency requirement of 1 sample. So no problem, we just read 1 sample in, did some processing within the sampling rate, and then read 1 sample out.

Now that this requirement has been relaxed, I think that we can use block processing to do things such as spectrum computation and processing and use more than 1 sample to do the adding and multiplication, etc.
Does this buy us blockSize * (1/fs) compute time?

My problem is, being a noob, I don't know how to formulate this new requirement in our chip's main() function.

So let's pretend we have some function SPI_Read(&mic) which is reading a sample, and SPI_Write(&output) which sends to a speaker. How would I transform this into a block based system?

The current system is setup like:

main() {
    while() { // check for interrupt
        // while loop executing on system timer set to 1 / fs
        SPI_Read(&in); // read from adc
        // filters
        SPI_Write(&out); // output to dac

    }
}

This is running on a cortex-M4, no codec.

Thank you for your help! (:


texts: Plotting the frequency response of a FIR filter

texts:
I am trying to figure out the mechanics of plotting the frequency response of a FIR filter. For example, I have used an algorithm (Parks-McClellan), to generate a low pass FIR filter with an even number of coefficients. I've got 34 tap values. Now I want to show the frequency response of this configuration in a normalized frequency 0 ... 1.0 (Fs) graph. The literature says, "The frequency response of the filter is computed by passing the array of coefficients through the discrete Fourier transform (DFT)." The text goes on to show this nice smooth graph, plotting magnitude vs frequency from 0.0 to 1.0. I would like to reproduce that graph.

So I sent it through my DFT code, well a 34 point DFT generates a very blocky DFT and it looks nothing like the graphic. I wanted more points of the DFT (between 0  and 34) and tried adding fractional points 0.0, 0.1, 0.2, .. 34.0) and got a very "loopy" DFT with lots of humps (and also not at all like the graphic), I tried "upsampling" adding zeros between the values, still no joy. I tried allowing the tap sequence to repeat, no joy. So I'm really curious how its done. Sending it to the MATLAB plot function produces the expected graph, so I have some confidence the data is good. Just my understanding of the mechanics of how to get the higher resolution graph are incomplete.


texts: Frequency response of numerical derivative

texts:
Analytical derivative of a function is equivalent to convolution of that function with $s$ in Laplace domain. Numerical derivatives are limited in bandwidth due to finite sampling rate, so they are not synonymous with convolving the signal with with $s$ term. At higher frequencies one would expect attenuation of the numerically differentiated signal from one that was computed analytically. Recently, I found that there are some differences at the low frequency limit as well which I cannot explain.

Attached is a plot of a signal sampled from a normal distribution (blue) and it's first derivative in time (red). As expected, at high frequencies the derivative signal begins to attenuate. But why does it not cross $\omega$ = 1 rad/s or 0.16 Hz as would be the case if the solution was obtained analytically?

Here's the code I am running in MATLAB

sr = 100000;
y = randn(1,sr);
dydt = y;
for i = 2:length(y)-1
    dydt(i) = (y(i+1)-y(i-1))*sr*2;
end
hold on, plot(abs(fft(y)));
plot(abs(fft(dydt)));
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')


texts: $\sin (t \omega)$ is not an Energy Signal, then how come its Fourier transform do exist?

texts:
The following integral (perhaps fourier tranform of $\sin (t \omega)$ ) is not convergent:

$\int_{-\infty }^{\infty } e^{-i t \omega } \sin (t \omega ) \, dt$

As, $\sin (t \omega)$ is NOT an Energy Signal (but a Power Signal), then how come we get successful in finding the fourier transform of $\sin (t \omega)$ ?


texts: Can I assume system is LTI when given by DTFT of impulse response

texts:
I'm having hard time to grasp it probably because i don't fully understand it.
I understand that when a system is given by $h(t)$ (in general $h(t-\tau)$) i can assume that it is a LTI system.
So i guess my gap of knowledge is that im not sure if it is possible to FT or DTFT impulse response of a non LTI system $h(t;\tau)$. I assume not because the transform depends on the time variable but here i have a 2 variables function and im not really sure how the transforms behave under these conditions.
So like the title says, if a system is given by $H(e^{j\theta})$ can i assume that it is a LTI system? and why is that?
Thanks in advance.


texts: the best pre emphasis for MMSE equalizer

texts:
The system consists of multiple low-pass filter and noise, $h_1$, $h_2$, $h_3$, ..., $n_1$, $n_2$, $n_3$, ...
At the transmitter, the power of signal is fixed; and at the receiver, MMSE equalizer is used to minimize the MSE.
My question is that if we can pre-emphasis the signal in the transmitter to improve the overall performance. If the distribution of signal and noise is known, what is the optimal pre-emphasis system?


texts: Interpolation FIR filter Remez design criteria

texts:
I want to upsample an audio signal by a factor of 8.

The technique I was going to use was to upsample by zero stuffing and then interopolate using a low pass anti-aliasing FIR filter designed using the Remez exchange algorithm.

My question is: what are the design constraints for the Remez algorithm?

My initial thoughts are:

If I have an initial sampling rate of 44100Hz then my new sampling rate will be 352800Hz. When I zero stuff, my spectral content in the range 0-22050Hz will be folder around the new upsampled Nyquist of 176400Hz in the range 154350-176400Hz. So I have the following constraints:

Sample rate: 352800Hz
Pass band: 0-22050Hz (0dB)
Stop band: 154350-176400Hz (-120dB)

This does not produce an very good filter, so my assumption that I only have one spectral image folded around 176400Hz must be incorrect.


texts: Nyquist criterion for zero ISI in the frequency domain

texts:
Can someone explain in simple words /graphs the necessary and sufficient Nyquist criterion for zero ISI in frequency domain, namely:

The constant folded Fourier Transform.


texts: Stability of open-loop transfer function from its Nyquist plot

texts:
I am facing a confusion on understanding system open-loop transfer function stability from its Nyquist plot.

According to the formula, for open loop transfer function stability: $$Z=N+P=0$$
where $N$ is the number of encirclements of $(0,0)$ by the Nyquist plot in clockwise direction.

Now how can I verify this formula for the open-loop transfer function: $$H(s)=-\frac{1}{s^3(s+1)}$$
The Nyquist plot is attached in the image.

Here $N=1$. Now how can I find $Z$ and $P$ and verify the formula?


texts:
I am doing a project on ECG arrythmia analysis using matlab.

I have designed notch filter for removing 50 Hz noise but don't know how to add a 50 Hz powerline interference noise to a clean ECG signal?
Also, I want to check whether noise is reduced in the filtered signal. Will Power spectral density using modified welch periodogram indicate whether noise is filtered or not?
How can I compare which wavelet (e.g. db6) is best suited for ECG analysis?


texts: Can somebody help me solve this signal accumulation problem?

texts:
Investigation of the method of suppression of random noise by coherent signal accumulation

Purpose - to identify opportunities for the coherent accumulation for cases of stationary and quasi-stationary signal.
Suppose that the input is the mixture of observed signal and a random white noise (i.e. noise with uniform spectral density). The signal is stationary and is described from sample to sample by constant function (for example, a sinusoidal signal with a constant frequency and constant initial phase). In this case, the input noise in it amplitude is several times greater than the amplitude of the signal. By coherent accumulation of the input mix for a number of samples it possible to increase the signal/noise ratio.

According to the results of modeling build a relationship:
a) The signal/noise ratio in the output mix with respect to duration of accumulation, ie accumulated number of samples at a constant signal/noise ratio at the input , ( the number of samples of accumulation varies )
b) the signal/noise ratio at the output with respect to signal/noise ratio at the input for a fixed number of samples (M = 10, 25 , 50 ) (SNR at the input vary)
Repeat item "1" for the case of quasi-stationary signals.
As a useful signal you should set a rectangular pulse of constant duration, which offset from the origin varies from sample to sample by a linear law.
To develop a functional diagram of a device that filters the signal by accumulation.

Type of signal: The sum of two harmonic signals;
S/N ratio: 0.2;
The number of cycles of accumulation: up to 500;
Limits of change in the signal / noise ratio: 0.1-2.

THEORY
The method of accumulation is applicable if the useful signal during the reception time constant or a periodic function. The method is consist of many times repeated signal and summation of its individual realizations in the  receiver.
Let the transmission of the desired signal implemented in two levels.

In the interval Tx signal is constant.
In the observation interval Tx sample values of the received signal are accumulated.

y1=x+r1
y2=x+r2
...............
yn=x+rn

and these values are summed.

We introduce two assumptions:
1) counts of interference ri are independent of each other;
2) The obstacle is stationary (its characteristics do not depend on time)
and lets find (Px/Pr)out on the output of the accumulator, ie

*вых -output


texts: chirp signal - Frequency spectrum characteristics

texts:
I am confused with my own intuition of the frequency response of a chirp signal.

In this case, the chirp signal is a frequency modulated signal with initial frequency starting at 100 Hz and ending at 300 Hz for 2 seconds.

Since the signal contains equal amounts of frequencies from 100Hz
to 300 Hz, the magnitudes of the frequencies should be equal? However, in the plot below, the frequencies at 100Hz and 300 Hz
have higher magnitude than the middle frequencies e.g. at 200Hz?
Should the plot be just a straight line?
Also, why does it appear that the middle frequencies are denser?
Is it correct to assume that since there are more frequencies in the
middle e.g. 100ish to 200ish than the end frequencies of 100Hz and
300Hz, then this should explain the increased density?

What is wrong with my assumptions?


texts: Equalization and then Sharpening vs Sharpening and then Equalization

texts:
This was a question from an exam of digital image processing:

"Working with low contrast images, an elective method to enhance details and improve contrast is to apply in cascade an equalization of the histogram and a sharpening filter technique.
a) Does the order of application of the methods matter? Justify
b) Which method would you use first? Justify"

Could anyone help me? I can't find a clear answer in the textbooks.


texts: Cross correlation vs FFT for finding phase between 2 signals

texts:
I have been struggling with this issue for months now, trying to understand which one to pick: FFT or cross correlation.

I am given a sine wave (with phase, amplitude, frequency and time vector known), and I have to extract the phase shift between that sine wave and another signal. The new signal has the same frequency and the same time vector as the sine wave, however the amplitude and the phase are unknowns. What would you guys suggest? I am at a loss.

Thank you.


texts: Compressive Sensing Through MATLAB Codes

texts:
I am new to the topic of compressed sensing. I read a few papers about it by R.Baranuik, Y.Eldar, Terence Tao etc. All these papers basically provide the mathematical details behind it, i.e., Sparsity, RIP, L1 norm minimisation etc. However, can anyone provide the MATLAB codes that achieve compressive sensing?

Thanks in advance for any help.


texts:
I've read that phase vocoder is the STFT, but done in windows (i.e. using windowing).

Why is it called "phase vocoder"?


texts:
I was recently experimenting with UKF that compensates nonlinear effects of noises by forming appended state vector. The general way of dealing with such case in the literature is the following:

Form appended state vector and expanded state covariance matrix.
Do a Cholesky decomposition (or some other square root decomposition) and form sigma points.
Propagate sigma points through state equation and calculate predicted state and state covariance matrix.
Additionally propagate predicted sigma points through measurement equation as well and form estimated output vector and output covariance matrix
Do standard UKF update calculations

This works fine, but I wondered why couldn't we do double Cholesky decomposition like we do in "vanilla" UKF with additive noises where we can do another sigma point generation with propagated state covariance matrix. So I tried to do first Cholesky decomposition where state vector was appended with state noises and expanded covariance matrix, and then do a second Cholesky decomposition with predicted state and covariance matrix which were appended by measurement noise effects. To my surprise the results of such filter were very bad. I cannot seem to explain to myself why would that be the case. Could anyone shed some light on why would such procedure not work very well.
Just as a note, I did make sure to adjust filter parameters taking into account sizes of appended state vectors (which were different in cases for first and second decompositions).


texts: strange phenomenon in spectrogram of Windowed Fourier Transform [edited]

texts:
I am trying to understand how the Windowed Fourier transform works, so I used Matlab to create an artificial signal, and plotted the Windowed Fourier transform using a gaussian window.

$x(t) = \left\{\begin{matrix}
 \sin(2\pi100t)& \text{for }0\leq{t}<\frac{1}{4}\\
 \sin(2\pi50t)& \text{for }\frac{1}{4}\leq{t}<\frac{1}{2}\\
 \sin(2\pi25t)& \text{for }\frac{1}{2}\leq{t}<\frac{3}{4}\\
 \sin(2\pi10t)& \text{for }\frac{3}{4}\leq{t}\leq1
\end{matrix}\right.$

and $f_s=500$, i.e. the data points are sampled at every 1/500 second.

This is what I have obtained for my plot.

However, for the last 2 intervals, $\frac{1}{2}\leq{t}<\frac{3}{4}$, and $\frac{3}{4}\leq{t}\leq1$, we can see some form of bifurcating phenomenon (vertically). Why is this so?


texts:
I am trying to formulate an n-dimension algorithm for the Cooley-Tukey FFT. In two dimensions, it is very clear:

$$
Y[k_1 + k_2 n_1] = \sum_{l_2 = 0}^{n_2} \Biggl[ \Biggl( \sum_{l_1 = 0}^{n_1 - 1} X \left[l_1 n_2 + l_2 \right] \omega_{n_1}^{l_2 k_1}\Biggl)\omega_n^{l_2 k_1}\Biggl]\omega_{n_2}^{l_2 k_2}
$$

Here is a link on more documentation on this formula: Implementing FFTs in Practice

In the case of 3 dimensions, how would this formula change? So if in the previous example we had $$ N = n_1 * n_2$$

How could I write the formula for $$ N = n_1 * ( n_2 *n_3 )$$

and in the case of n dimensions, how could I general the summations? I am having trouble understanding how the twiddle factors get accounted for.


texts:
Assume I have a square wave pulse train, where each high lasts for H and each low lasts for L seconds. I sample this waveform and have the sequence of numbers.

Question is what should be the sampling rate if I just want to know what level I am in(High or Low)?

I know by Nyquist I can find the sampling rate for full-reconstruction of the pulse train, but I don't need the full reconstruction. I just want to be able to:

1. Don't miss any High or Low

2. Detect all Lows and Highs


texts: Some difficult-to-approximate signals to test an signal extracting algorithm

texts:
I am searching for some strange signals in order to test an algorithm that I am developing. I've found an article by Donoho and Johnstone where they propose to test the algorithms using Blocks, Bumps, Heavysine, and Doppler signals.
Is there any other signal that is generally difficult to be estimated? I am not really into this topic and I would really appreciate some help.
Thanks in advance.

Edit: As suggested in the comments I'll provide some context. I am working with Spline regression in noisy context. My idea is to provide an alternative solution for the Smoothing Spline by using an adaptive algorithm that chooses the location of the knots according to the behaviour of the data. For instance if I would like to approximate a Doppler signal I would expect the algorithm to put very dense knots on the lefter side to capture the rapid movements and sparse knots in the righter side instead.
I've been suggested to ask this in here even though I am not quite sure if it's a suitable question for this particular site.


texts: Simulation of images of a fluorescent bead

texts: Figure 2 in the paper by Thompson et al. "Precise Nanometer Localization Analysis for Individual Fluorescent Probes" (https://pubmed.ncbi.nlm.nih.gov/11964263/) presents series of images of a fluorescent bead for different numbers of photons N in the spot. Could someone give me some practical hints how to simulate such images? I see that the resolution of each image is the same and the quality depends on the numbers of photons N. The authors explain that "images with different numbers of photons were created by accumulating, offline, varying numbers of frames from a time-series stack of a single bead" but it is still not clear for me how to proceed.


texts: Is it possible to specify number of taps and cutoff frequency independently?

texts:
When specifying the properties of a moving average filter, is it possible to specify the number of taps and the cutoff frequency?
My guess is no since taps (the number of samples in a moving average filter) is what determines the cutoff frequency according to the frequency response equation.
If I'm wrong, how is it possible to specify the two?


texts: KLT tracker, reject tracked occluded points

texts:
I want to use KLT tracker for a visual odometry application. Thus, I only want to track object features. Currently I'm using OpenCV's implementation  cv.calcOpticalFlowPyrLK().

I'm extracting Shi-Tomassi goodFeatures2Track and doing bi-directional checking (i.e. points tracked from frame a to b must be tracked from b to a as well), which filters several outliers.

However, when two objects intersect due to occlusion or reflections, a point is detected and tracked. these points usually slide along the edges of objects as frames pass, this makes sense given that KLT minimizes an appearance transform. However, they don't describe the structure of the environment being tracked.

I'm able to filter some of them using 5-point epipolar RANSAC for tracked points with long enough displacement, but tome of them remain.

Does anybody know how to filter these sliding tracked outliers?


texts: CDF of product of two *translated* exponential variables

texts:
It is known that the CDF of the product of two exponential random variables is given by this answer

  If $X$ and $Y$ have rate parameters $\mu$ and $\lambda$, so their
  pdf's are $f_X(x) = \mu e^{-\mu x}$ and $f_Y(y) = \lambda e^{-\lambda
> y}$ for $x,y \ge 0$, then for $c > 0$, $$ Pr(X Y < c) = \int_0^\infty
dx \int_0^{c/x} dy\; \mu \lambda e^{-\mu x - \lambda y} = 1-2\,\sqrt {
{c\mu\lambda}}{{ K}_{1}\left(2\,\sqrt { c\mu\lambda}\right)} \tag{1}$$  where
  $K_1$ is a modified Bessel function of the second kind.

I need to calculate the CDF of $Z = UV$ where $U=a+X$ and $V=b+Y$ that $a,b>0$ are constant.

It is easy to show that the PDF of $U,V$ are $f_U(u)=\mu e^{a\mu} e^{-\mu u}$ and $f_V(v)=\lambda e^{b\lambda} e^{-\lambda v}$ for $U \geq a > 0$, $V \geq b > 0$.

$$ Pr(Z < c) = \int_a^\infty du \int_b^{c/u} dv\; \mu \lambda e^{a \mu + b\lambda} e^{-\mu u - \lambda v} \tag{2}$$

I don't know how to process further.


texts: How do wavelet levels depend on one another?

texts:
I am working with wavelets for my thesis, and I would like to know if there is an intuitive dependency between the levels of a DWT.

To make it clear: The DWT performs a convolution, and then the result is downsampled, right?

Like in this example of a convolution below.

Now when we downsample dyadically, we throw away every second value, right?
Like this:

Can we say that 2 adjacent values in the original are projected onto 1 value in the transform?

For example, the first value in the first level of the DWT would then be "responsible" for values 1 and 2 in the original.
So that we would know where e.g. a pattern in a high DWT level actually comes from.

Does this kind of projection exist?
If it does not work, why not?
Is the convolution/downsampling performed differently than I
imagined?


texts: Why is the median of a PDF not always equal to the mean?

texts:
When studying probability density functions(PDFs), the focus is often placed on the "mean" and not the "median". The mean is often described as the center of mass of the PDF. I've always taken this as true.

However, the definition of the median seems to correspond to what I would consider the center of mass. That definition being:

$\int_{-\infty}^{median}pdf(x)dx = \frac{1}{2}$

Looking at this, we're essentially integrating the PDF until we reach the half-way point. To me, this seems identical to what we would consider the center of mass.

So why does integrating a PDF until we reach halfway (the median) not always equal the mean if the mean is the center of mass?

Note that I understand the difference between the mean and median from a general statistics perspective, but I'm curious from a PDF perspective.

I'll also format the question in the following manner:

$\int_{-\infty}^{\infty}pdf(x)dx = mass$

$\int_{-\infty}^{\infty}x*pdf(x)dx = center\:of\:mass$

If we define the following equation:

$\int_{-\infty}^{T}pdf(x)dx = \frac{mass}{2}$

Why is T not always the center of mass? It appears T is the center of mass if the function is symmetric, but why doesn't it hold true if the function is non-symmetric?


texts: MU-MIMO Sensitivity to Line of Sight (LOS)

texts:
Why do researchers talk about MU-MIMO being "less sensitive to line of sight" / "LOS Okay"? For example in this video, time=[00:40,00:42]

Similarly, "User separation based on spatial channel properties is
particularly  difficult  in  situations  where  the  users  are  located
close to each other and experience LOS propagation conditions
to the BS antenna array"

My perception is that LOS is always preferred to the opposite, multipath propagation.

Can someone please explain this? Thank you.


texts: Demodulation distortion in single-sideband and double-sideband

texts:
I am trying to compare the demodulation for single sideband (SSB) and dual sideband modulation formats (DSB) when the reconstructed carrier is not perfect (i.e. shifted in both frequency and phase). I want to show mathematically in which case distortion would be more problematic.

For instance, if the carrier is $\cos (2 \pi \nu_c t)$, the shifted reconstructed carrier would be:

$$\cos (2 \pi (\nu_c + \Delta \nu)t + \varphi) \tag{1}$$

Attempt:

SSB: A signal transmitted with the single sideband (SSB) modulation is given by:

$$f(t)=A\ \Re\left\{m_a (t) e^{j2\pi \nu_c t}\right\} = A\left[m(t) \cos(2\pi \nu_c t)- \hat{m}(t) \sin (2\pi \nu_ct)\right] \tag{2}$$

where $m_a$ is the analytic signal and $\hat{m}(t)$ is the Hilbert transform. Demodulation is obtained by multiplying the transmitted signal with the reconstructed carrier given in $(1)$. With some manipulation this became:

\begin{align}
A\left[\frac{m(t)}{2}\left[\cos [(2\omega_c + \Delta \omega)t+ \varphi]+ \cos(\Delta \omega t + \varphi)\right]\\
- \frac{\hat{m}(t)}{2}[\sin[(2\omega_c+\Delta \omega)t+\varphi]+\sin(\Delta\omega t +\varphi)]\right]
\end{align}

Low pass filtering:

$$\therefore \frac{A}{2}\left[m(t)\cos(\Delta \omega t + \varphi)- \hat{m}(t)\sin(\Delta\omega t +\varphi)\right] \tag{i}$$

DSB: And the dual sideband modulation (DSB) signal is given by

$$f(t)=Am(t) \cos(2\pi \nu_c t) \tag{3}$$

Demodulation gives:

$$Am(t) \cos(2\pi \nu_c t) \cos (2 \pi (\nu_c + \Delta \nu)t + \varphi)=\frac{A m(t)}{2}\left[\cos((2\omega_c + \Delta \omega)t + \varphi)- \cos (\Delta \omega t + \varphi)\right]$$

Low-pass filtering:

$$\therefore -\frac{A m(t)}{2} \cos (\Delta \omega t + \varphi) \tag{ii}$$

So, if my approach is correct, how can we interpret equations (i) and (ii) to see which modulation format results in having more distortion?

Any explanation would be greatly appreciated.


texts: IFFT of a complex, conjugate-symmetric transfer function - Value at Nyquist frequency

texts: Using the MATLAB 'symmetric' flag to compute the IFFT of a complex-valued, length-N (with N even), transfer function, the imaginary part of the complex value at the index ((N/2)+1) corresponding to the Nyquist frequency appears to be ignored. Why is this done instead of taking the magnitude of the complex value at the Nyquist frequency? How important are the values at the Nyquist frequency and DC for audio applications?


texts:
Is there a simple way to compensate for the reduction in amplitude of a low pass filter of a noisy sine wave signal that has been smoothed using a FIR filter so that the amplitude of the smoothed sine wave matches that of the approximately known amplitude of the noisy original signal? The approximate amplitude of the original data could be estimated empirically from the RMS value over one full (estimated) period. The FIR filter coefficients I am using are

h = [ 1 2 3 3 2 1 ] / 12


texts: How to convert filter into frequency domain to do filtering?

texts:
I have this filter:

filter_2 = firceqrip(2,0.6,[0.05 0.03]);

I want to convert it to the frequency domain to multiply it by a signal (i.e filter in the frequency domain)

I have a signal and I have a filter.

I converted the signal into frequency domain and now I need to filter it with my filter.

Converting the signal into frequency domain is easy, but how do I filter the signal now?

This is my signal:

freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
t_max = 4 * periods(1);
t = linspace(0, t_max, 50);

signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);

% to find the fft of the signal
N = 64;
signal_spect = abs (fft(signal,N));
signal_spect = fftshift(signal_spect);
F = [-N/2:N/2-1]/N;

now how do I filter the previous signal with my filter?
I tried doing this fft(filter_2) but it didn't work out, it's probably a wrong approach, but I do not know what else to do

any help is appreicated


texts: How do I properly do gamut mapping for an image?

texts:
I asked this question earlier. To summarize, dithering an image with a specific palette was causing the image to look smeared and very discolored.

Thanks to Olli Niemitalo, I found out my issue was because the colors in the original image (particularly the near-blacks and near-whites) were outside the gamut of the palette I had. As a result, Floyd-Steinberg dithering was piling up tremendous error levels which were being carried over to the next pixels. I confirmed this to be true by replacing every out-of-gamut color with the nearest color from the palette (which are guaranteed to be in the gamut, by definition). Trying to dither after doing this yielded better results.

What I'm having trouble with is finding the nearest color inside the gamut of my palette, given another color. Here's my current method:

Convert palette colors from sRGB to XYZ.
Interpret them as a collection of 2D points on a grid (discarding Z)
Draw a convex polygon around the points to determine the gamut.
For each pixel in the image

Convert the pixel from sRGB to XYZ
If the pixel's XY values are not in the polygon, find the nearest point that is inside it.
Convert this result back to sRGB (using the new X and Y values, and the original Z value)
Replace that pixel in the image with the new color

After all that, dither with Floyd-Steinberg.

This is the palette I'm working with (in (R,G,B) format):

(162, 143, 66)
(128, 148, 168)
(120, 99, 100)
(31, 39, 97)
(126, 116, 103)
(203, 35, 9)
(57, 81, 43)
(104, 101, 98)

Original image:

Gamut mapping (attempt):

Dithering:

That certainly looks better than the original problem I was having, but it's still troublesome. I'm sure I'm getting the gamut mapping procedure wrong somehow, and I suspect it has to do with discarding Z, but I have no idea how to proceed. Googling this yielded only Photoshop instructions, which isn't what I'm looking for.


texts: Which domain used for denoising additive and multiplicative noises

texts:
I want to know what is the difference between additive noise, multiplicative noise.. In what domain these noises are handled to remove?

I am interested particularly in image dataset

Does different domains are used in case of satellite imaging where noise due to reflection of frequencies is statistically independent?


texts: I am having this simple function which acts on two state variable $s_t$ (states being $\{0,1\}$, boolean in sense), and I am evaluating it in discrete time steps. The function is $f_t=1$ if $s_t$ has been $1$ for last $W$ samples (some fixed window size), and $f_t=0$ otherwise. So basically this will set the output to active if the input signal has been active for certain amount of time/samples. Is there some name or category for this? I thought it is latching but according to definitions on web it seems latching is more about delaying the signal, instead of asking "was the signal delayed on input".


texts: Upsampling - What purpose does the interpolation filter have?

texts:
I want to apply some nonlinear processing to a signal, namely: I want to implement a tube emulation which adds warmth/harmonic distortion to a digital audio signal. I am worried about aliasing, so I figured out the following processing sequence:

Upsampling of the sample rate from $f_s$ to $2 \cdot f_s$
Applying an interpolation lowpass filter with $f_g = 0.5 f_s$ and a passband gain of $2$
Applying the nonlinear processing to the upsampled signal
Again applying a filter with with $f_g = 0.5 f_s$
Decimation back to $fs$.

My question arises at point 2, the interpolation filter. If I google this topic, it results in fascinating concepts like Lagrange interpolation, spline interpolation and polynominal interpolation and the math behind those seem far from trivial. Moreover, I don't understand why this is even needed - what is the reason why one would utilize one of those concepts rather than just running the signal through a run-off-the-mill IIR with sufficient stopband attenuation? The only thing I can understand is why the passband gain needs to be non-unity.

Also, although not main point of my question, feel free to point out any logical flaws that might be in my processing sequence. Thank you!


texts: Relationship between z-transform and DFT

texts:
I'm studying for a Signals Processing exam and came across an exercise that I'm finding pretty difficult to solve. It says:

  Asume there is a signal $x[n]$ of length $N$. Its $\mathcal{Z}$-Transform is $X(z)$. Find the expressions of the DFTs of the next signals using a
  correct sampling of $X(z)$.

  a) $x_{1}[n]=x[n]$  if  $0\le n \le N-1$ and $x_{1}[n]=0$  if  $N\le n \le 2N-1$

  b) $x_{1}[n]=x[n]+x[n-N]$

For a), I thought of expressing  $X[k] = X(z) \bigg|_{z=e^{j\frac{2\pi}{N} k}}$, i.e. the $\mathcal{Z}$-Transform evaluated in $z=j\frac{2\pi}{N} k$ is the DFT of $x[n]$. Due to the zeros added at the end of the original signal, I arrived at the conclusion that  $X[k] = X(z) \bigg|_{z=e^{j\frac{2\pi}{2N} k}}$. Namely, the $\mathcal{Z}$-Transform is now sampled in $2N$ points around the unity circle. Is this reasoning right?

And regarding exercise b), I have absolutely no idea how to solve it.

Thanks for your time!


texts: Why does adding a time delayed version of a signal to itself create a filtered signal?

texts:
I was asked this question and couldn't come up with an answer on the spot which didn't involve the frequency domain (basically that the co-efficients of the delay sequence are the impulse response of a FIR filter).

Does anyone have any insight which makes this process 'obvious'?


texts: How to explain this BER vs SNR behaviour?

texts:
I wrote a script in MATLAB to simulate an LDPC created using an algorithm of a colleague of mine. The parity check-matrix is regular and has dimensions $1012 \times 1518$ and the code rate is $R=1/3$.

The channel I'm using is a BPSK one with added Gaussian white noise. The mapping of the bits I used were the values $-1$ and $1$. The decoding method is sum-product (or belief propagation).

To perform the simulation, I used as a guide the Example 2.6 in the 36th page of this paper. I'm getting this curve:

This is really different from what I would have expected. The fact that the function is no monotonically decreasing was very shocking. Also, the fact that there is a BER of approximately $1$ for a positive SNR doesn't make any sense to me.

I'm thinking that, maybe, what is wrong is how I defined the a priori LLRs or the SNR. Given that I used $E_b = 1$ to map the bits, then I used some formulas from the section 3.2 of this paper:

$$R\cdot E_b = E_s$$
$$N_0 = 2\sigma^2$$

So, using the formula mentioned in the Example 2.6 of the first paper, we get that

$$r_i = 4y_iR\frac{E_b}{N_0}$$

I believe that this weird behaviour is coming from that formula, maybe a difference of definitions between the two papers that I didn't take into account. I don't really know, what do you think is the reason of the function having that shape?

Extra info:

The SNR axis in the plot corresponds to $R\frac{E_b}{N_0}.$
I've read here a question that could solve my problem, but I'm not sure if it applies directly. Here, I don't have any information about the data bit rate that is mentioned in that question. I don't really know what that quantity they call $R_b$ stands for, either.


texts: Modifying sizes and thickness of letters to "embed" an image in text

texts:
I am wondering if this is defined, names and possibly solved - suppose I have a grayscale image and about one page worth of text. I would like to attempt to give  the text an appearance of an image when looked from far away. This would mean that text would be mapped to the space of an image and each letter would have size and thickness proportional to the mean pixel value of the area it maps onto. Any output  is ok (vector or raster).

Is this something that is used, named, does anyone know a place to look for this before I start writing the code from scratch? Python, R, MATLAB or a command line tool would all be acceptable.

If the whole pipeline is not a known use case, I would appreciate a relatively easy way to put letter X with (locX,locY,size,thickness) = (x,y,s,t) in an image, and I would write the rest myself.

(This question is too general for stackoverflow, but if it should be moved elsewhere I would appreciate that help as well.)


texts:
I have two Time Domain functions, $f_1(t)$ and $f_2(t)$.

I have both Fourier Transforms,  $F_1(\omega)$ and $F_2(\omega)$. Functions $f_1$ and $f_2$ are not independent and, in fact, $f_1$ is also a function of $f_2$. So, the derivative $f_3(t)=\frac{d\,f_1}{df_2\,}$ is meaningful.

Can I obtain the Fourier Transform $F_3(\omega)$ (the FT of $f_3$) directly from $F_1(\omega)$ and  $F_2(\omega)$?

I wish to calculate it directly in the Frequency Domain, without going to the Time Domain to do the derivative.