texts: What is the best way to show beamformer enhances signal of interest over interference

texts:
I need to study an array's performance in terms of interference.
I have already analyzed the array's white noise gain which is the ratio of output and input SNR of the beamformer i.e. $\rm{SNR_o} / {SNR_i}$. This was based on defining array's output as
$${\bf y}(t, \phi) = {\bf a}(\phi)s(t) + {\bf n}(t),$$
where $s(t)$ is the signal of interest, ${\bf a}(\phi)$ is the array manifold, and ${\bf n}(t)$ an isotropic AWGN that is spatio-temporally uncorrelated (and uncorrelated with $s(t)$). Note that ${\bf y}(t, \phi), {\bf a}(\phi), {\bf n}(t) \in \mathbb{R}^{M}$, $M$ being the number of sensors in the array.

If an interference signal is added, the array outputs
$${\bf y}(t, \phi_s, \phi_v) = {\bf A}(\phi_s, \phi_v){\bf r}(t) + {\bf n}(t),$$
where ${\bf A}(\phi_s, \phi_v) := [{\bf a}(\phi_s) \;  {\bf a}(\phi_v)]$,
${\bf r}(t) := [s(t)\;  v(t)]^T$, $s(t)$ is the signal-of-interest and $v(t)$ is the interference signal. $\phi_s$ and $\phi_v$ are the angles of arrival of the signal of interest and interference, respectively.

The idea is to study how the beamformer improves the signal of interest over the interference.

The concept of signal-to-noise-plus-interference ratio ${\rm SNIR}$ is introduced.

Question 1: Is the ${\rm SNIR}$ the best choice of metric for this study?

At the beamformer's input,
${\rm SNIR}_i = P_s / (P_v + P_n)$, ${\rm SNR}_i = P_s / P_n$ and ${\rm SIR}_i = P_s / P_v$. The interference-to-noise ratio ${\rm INR}_i = P_v / P_n$, where $P_s, P_v, P_n$  are signal-of-interest power, interference signal power, and additive noise power, respectively.

${\rm SNIR}_i = \dfrac{P_s}{P_v (1 + {\rm INR}_i^{-1})} = \dfrac{{\rm SIR}_i}{1 + {\rm INR}_i^{-1}}$.

The beamformer weight is ${\bf w}$, and the beamformer outputs
\begin{eqnarray}
B &=&
{\bf w}(\phi_L)^T ~ {\bf y}(t, \phi) \nonumber \\
%&=&
%{\bf w}(\phi)^T {\bf A}(\phi){\bf x} \nonumber \\
%&=&
%{\bf w}(\phi)^T \left[ {\bf a}(\phi_s)s(t) + {\bf a}(\phi_v)v(t)  + {\bf n}(t)
%\right] \nonumber \\
&=&
{\bf w}(\phi_L)^T{\bf a}(\phi_s)s(t)  +  {\bf w}(\phi_L)^T {\bf a}(\phi_v)v(t)  +{\bf w}(\phi)^T{\bf n}(t) \nonumber
\end{eqnarray}
where $\phi_L$ is the beamformer's look direction.

Assuming the signal of interest, interference, and noise are uncorrelated,
the output signal-to-noise-plus-interference ratio
\begin{eqnarray}
{\rm SNIR}_o &=&
\dfrac{
 {\bf w}(\phi_L)^T {\bf a}(\phi_s) {\bf a}(\phi_s)^T {\bf w}(\phi_L) ~P_s
}{
 {\bf w}(\phi_L)^T {\bf a}(\phi_v) {\bf a}(\phi_v)^T {\bf w}(\phi_L)~ P_v
+
 {\bf w}(\phi_L)^T {\bf w}(\phi_L) ~P_n
} \nonumber \\
&=&
\dfrac{
  {\bf w}(\phi_L)^T {\bf a}(\phi_s) {\bf a_s}(\phi)^T {\bf w}(\phi_L) ~ {\rm SIR}_i
}{
  {\bf w}(\phi_L)^T {\bf a}(\phi_v) {\bf a}(\phi_v)^T {\bf w}(\phi_L)
 +
  {\bf w}(\phi_L)^T {\bf w}(\phi_L) ~{\rm INR}^{-1}
}
\end{eqnarray}

And the beamformer's ${\rm SNIR}$ gain
\begin{eqnarray}
\dfrac{{\rm SNIR}_o}{{\rm SNIR}_i} &=&
\dfrac{
  {\bf w}(\phi_L)^T {\bf a}(\phi_s) {\bf a}(\phi_s)^T {\bf w}(\phi_L) \left[1 + {\rm INR}^{-1}\right]
}{
  {\bf w}(\phi_L)^T {\bf a}(\phi_v) {\bf a}(\phi_v)^T {\bf w}(\phi_L)
 +
  {\bf w}(\phi_L)^T {\bf w}(\phi_L) ~{\rm INR}^{-1}
}
\end{eqnarray}

For the beamformer, $\underset{\phi}{\arg\max} ~{\bf w}(\phi_L)^T {\bf a}(\phi) = \phi_L$, where $||{\bf w}|| = 1$. Assume that $\phi_s = \phi_L$ while $\phi_v = \phi_L \pm \Delta_{\phi}$.

Question 2:
At this point, what is the best way to show that this beamformer enhances the signal-of-interest over the interference signal?

Question 3:
Is it acceptable to do this study using only signal-of-interest and interference while ignoring the additive noise? That way, I will have just ${\rm SIR}$ gain instead of ${\rm SNIR}$ gain to consider.


texts: I am working on a project that involves time synchronization and we want to derive worst case bounds for the accuracy of the current time given measurements with known error as input. Unfortunately all the PLL literature I have found focuses on small-signal cases and transient response rather then error propagation. I see some reference to more statistically grounded techniques but haven't really found a great source. Any ideas where I might look next?


texts: Decimation function implemented in standard C

texts:
I'm thinking in a decimation C function. The prototype of the function is:

decimate(int N, float input[], float r, float output[])

where N is the length of array input and r is the reduction factor of sampling rate.

Is there a simpler way that decimation (matlab function) to implement a decimation fuction?

Somebody knows an Open Source C library for this porpouse?

Thanks.


texts: Symbol synchronization with M-QAM. M < 16

texts: I working on M-QAM system using a Cortex-M4 processor. Using a alternating symbol as a preamble, I am able to get carrier lock. Symbol synchronization is then achieved by finding the peak. An early late gate is used to adjust for clock differences between modulator and demodulator. An integrator is used to even out the noise from the error detector. However, I can see the peak slipping from the optimum sampling position. Eventually after 500-600 symbols it loses sync. Are there any superior ways of doing symbol sync? Am I missing something?


texts:
I am looking for an efficient head pose estimation algorithm from facial features.

What I already have: I already have implemented the code for stable facial landmarks(68 landmarks on the face).
I also have the intrinsic camera parameters.

The task is to efficiently determine the head pose from webcam video i.e,
(roll, yaw and pitch).

Which are the viable algorithms for the purpose?

I am working on python by the way, any working code for the same would be welcome.


texts: SPL dB to voltage conversion misunderstanding

texts:
base on these links:

https://electronics.stackexchange.com/questions/96205/how-to-convert-volts-to-db-spl
https://electronics.stackexchange.com/questions/302833/convert-digital-and-analog-values-to-db-spl
https://electronics.stackexchange.com/questions/354313/volts-to-decibels-using-microphones-transfer-function

they all say that sensitivity is reported in RMS but this link suggests:

https://www.analog.com/en/analog-dialogue/articles/understanding-microphone-sensitivity.html#

One last very confusing piece of this comparison is the inconsistent usage of peak and rms levels between digital and analog microphones. The microphone’s acoustic input levels in dB SPL are always rms measurements, regardless of the type of microphone. The output of analog microphones is referenced to 1 V rms, as rms measurements are more commonly used for comparing analog audio signal levels. However, the sensitivity and output level of digital microphones are given as peak levels because they are referred to the full-scale digital word, which is a peak value. In general, this convention of using peak levels to specify the output of digital microphones must be kept in mind when configuring downstream signal processing that may rely on precise signal levels. For example, dynamic range processors (compressors, limiters, and noise gates) typically set thresholds based on rms signal levels, so a digital microphone’s output must be scaled from peak to rms by lowering the dBFS value. For a sinusoidal input, the rms level is 3 dB (the logarithmic measure of (FS√2) below the peak level; this difference between rms and peak may be different for more complex signals. For example, the ADMP421, a MEMS microphone with pulse-density-modulated (PDM) digital output, has a sensitivity of –26 dBFS. A 94 dB SPL sinusoidal input signal will give a –26 dBFS peak output level, or a –29 dBFS rms level.

that sensitivity is reported in peak voltage not RMS.
which one is correct?


texts:
I am not sure if it's DSP question or physic.

But I wonder is there any math equation (method) that describe that graph? (It could be anything approximately like that).

I am working on audio freq analyser. And for typical pop song I get graph that Low End is very big but Middle and High End are very low. That's the mathematical reality, but I want to transform it to human ear reality. Could anyone help me?

(Of course I know how to make freq log scale on horizontal axis and dB scale on vertical axis, just have problem to modify it to human ear characteristic)


texts: Two-GSMK channels demodulation: 3dB difference in BER

texts:
I am trying to demodulate individually two GSMK channels separeted by 50kHz (one channel = 25kHz).
I am working on a limited ressources hardware so I decided to make some optimization on the signal processing chain.

Here is the architecture: after sampling, my two channels are aliased at 4.6MHz and at 4.65MHz.

The architecture for one channel is as follow: a DDS (Direct Digital Synthesis) and a mixer bring my signal in baseband. I then have a downconversion, filtering, the decoding algorithm and an error controler to compute a BER.

The basic architecture for demodulating two disctinct signals would be to duplicate this processing. Unfortunately, my hardware cannot contain so much computing.

One of the idea I had was to put in common some DSP operations. The downconversion uses a lot of ressources so I decided to mutualize the process for both channels.
The downconversion uses a CIC filter, it filters and decimates the signal.

The new signal path would be: DDS+MIXER (channel 1 brought to baseband, channel 2 centered around 50kHz), CIC + compensation filter (CIC response is not flat in passband), second dds + mixer to bring channel 2 to baseband.

here is a scheme of the processing:

The "tricky" thing is that on my second DDS I multiply both I and Q by cosine instead of cosine / sine (in that case, it does not work).
When doing this, I guessed that maybe I would lose something.

When measuring the BER, I can easily see that I have a 3dB loss between my two channels. My question: how can I "mathematically" explain this?

Is it because I am mixing two cosine instead of one cos and one sine; is half of the information lost?
Am I losing 3 dB on the SNR?

Edit: Here are the BER curves


texts:
There are tons of papers discussed this topics and most of them are related to fancy optimization techniques (convex/non-convex). I am wondering if we can have a simpler way to find the antenna weights if I restrict my requirement as following:

Linear Array with antenna positions known at hand
All antennas are located on a grid. For example, four antennas are located at (0, d, 4d, 6d) where d is a fixed distance (say, half wavelength)
For synthesized pattern, we only care about mainlobe width and sidelobe level.

Any references are welcome

Thanks


texts: Filter a negative SNR signal whose frequency spectrum lies inside the noise frequency spectrum

texts:
I am a newbie in Signal Processing and need your suggestion in the following problem. My signal amplitude is 10 times weaker than the noise. Furthermore, the frequency spectrum of the noise and the signal is such that the PSD of the signal falls inside the PSD of the noise window (so convolution or frequency domain filtering is not working at all).
In the figure below, I show the time domain evolution of the signal, where the solid black line is "Signal + Noise". The dotted black line is "Only Noise". The red line represents the "actual signal".
Can you please suggest me which filtering method should I try to know whether the signal is present or absent in the noise? I am trying to implement this in Python.

And the PSD looks like this (Blue = Noise, Orange = Signal)


texts: Average image over $N$ frames with one frame buffer

texts:
Can you give some filter that is averaging an image over $N$ samples in a recursive way that only needs one frame buffer?

$\frac{1}{N} \sum_{k=1}^{N} y_{k}(m,n)$

I can only imagine filters that need $N-1$ frame buffers.


texts:
Any way to reverse engineer professional audio algorithms or those in e.g. rack units?

In order to study how they're implemented?

Not suggesting to steal ideas, but as education.


texts:
I am really confused about oversampling and upsampling. I know that upsampling means stretching our signal in time domain by a factor of K, then interpolating some values or using just zeros. As I found out stretching in time domain is shrinking in frequency domain by the same rate. For example, stretch by factor of 2 in time domain mean shrinking by factor of 2 in frequency domain.

Let talk about oversampling, we just increase the sampling frequency more than Nyquist rate. Then we have too much information. Is there any chance to have aliasing? I mean in time domain, lets say I am sampling my signal (frequency is 70 Hz) by sampling frequency equals to 280 Hz (Nyquist*2) so, is there aliasing in frequency domain? Sorry if my explanation is not clear. It is because I am bewildered with these concepts; oversampling, upsampling, downsampling, undersampling, critical sampling and aliasing.

theoretically, I know them but when it comes to practical view, I have serious problem.


texts: What is the intuitive idea behind FIR filter order and filter parameters?

texts: I know that if we decrease stop-band and pass-band attenuation, and also decreasing the transition band, the filter order will increase as well as filter performance. But what is the intuitive idea behind that? Thank you very much.


texts:
How do I calculate the Spectral Entropy of a signal in MATLAB ? I know the basic steps but it would be nice if someone can help,

Calculate the power spectrum of the signal using FFT command in MATLAB.
Calculate the Power Spectral Density using the power spectrum or using any other technique.
Normalize the Power Spectral Density  between $[0, 1]$, so that it can be treated as a probability density function $p_i $.
Calculate the Entropy $H(s) = -\sum p_i\log_2\left(p_i\right)$


texts: Wrong estimation of derivatives with an extended Kalman filter

texts:
I am trying to implement an extended Kalman filter (EKF) in MATLAB for the estimation of joint trajectories (angular position, angular velocity and angular acceleration) from noisy motion capture measurements using a constant angular jerk process model. This is my process model:

$x_{prior}^k = \begin{bmatrix} 1 & \Delta t & \frac{\Delta t^2}{2} & \frac{\Delta t^3}{6} \\ 0 & 1 & \Delta t & \frac{\Delta t^2}{2} \\ 0 & 0 & 1 & \Delta t \\ 0 & 0 & 0 & 1 \end{bmatrix} \cdot x_{postiror}^{k-1}$

and process noise covariance:

$Q = 2 \cdot \sigma_p^2 \cdot \begin{bmatrix} \frac{\Delta t^7}{252} & \frac{\Delta t^6}{72} & \frac{\Delta t^5}{30} & \frac{\Delta t^4}{24} \\
\frac{\Delta t^6}{72} & \frac{\Delta t^5}{20} & \frac{\Delta t^4}{8} & \frac{\Delta t^3}{6} \\
\frac{\Delta t^5}{30} & \frac{\Delta t^4}{8} & \frac{\Delta t^3}{3} & \frac{\Delta t^2}{2} \\
\frac{\Delta t^4}{24} & \frac{\Delta t^3}{6} & \frac{\Delta t^2}{2} & \Delta t \end{bmatrix}$

My estimates for the angular positions look ok, but my estiamtes for the derivatives are quite wrong. Here is a graph from a simple example application with a comparisson between the angular velocity estimated by my implementation (EKF) and a quite noisy central difference (CD).

Somehow, I cannot improve the estimation of the derivatives by tuning the covariances, initial values or (almost) anything else. So I think there must be something wrong with my implementation. But I cannot find anything, maybe I suffer from tunnel vision...

I adopted the EKF equations from the paper by Yu et al. 2.

% Apply extended Kalman filter equations
for k = 1:T

% Apply time update
if k == 1

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(x0);
    A = dfdx(x0);
    W = dfdw(x0);

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P0 * A' + W * Q * W';

else

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(X_posterior((k - 1), :)');
    A = dfdx(X_posterior((k - 1), :)');
    W = dfdw(X_posterior((k - 1), :)');

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P_prior(:, :, (k - 1)) * A' + W * Q * W';

end

% Apply measurement update
H = dhdx(X_prior(k, :)');
V = dhdv(X_prior(k, :)');

% Equation (11) in [Yu2004] using the pseudo inverse to handle
% badly scaled matrices
K = P_prior(:, :, k) * H' / (V * R * V' + H * P_prior(:, :, k) * H');

% Equation (14) in [Yu2004]
X_posterior(k, :) = X_prior(k, :) + (K * (Z(k, :)' - h(X_prior(k, :)')))';

% Equation (10) in [Yu2004]
P_posterior(:, :, k) = P_prior(:, :, k) - K * H * P_prior(:, :, k);

end

There seems to be one difference in the computation of the prior error covariance compared to the equations given in the paper by Welch and Bishop 3.

Yu et al.:

$P_{prior}^k = A \cdot P_{prior}^{k - 1} \cdot A^T + W \cdot Q \cdot W^T$

Welch and Bishop:

$P_{prior}^k = A \cdot P_{posterior}^{k - 1} \cdot A^T + W \cdot Q \cdot W^T$

The equation by Yu et al. produces the wrong derivatives and the equation by Welch and Bishop seems to result in instable estimations.

I compiled a simple example application with my implementation 4. It should work out of the box, if you run "TestKalmanFilter.m" in MATLAB and you should see the given graph and a small animation of a double pendulum.

Any help is greatly appreciated.
Janis


texts:
I am trying to make a DFT on a signal with frequency f=50 MHz main component plus some noise.
As far as I know if I sample it at F=100MHz I should be able to get a proper plot of the DFT since F=2f and I should see a peak at 50MHz, this instead does not happen,see first figure.

But instead I obtain what expected when I am leaving unchanged all parameters of DFT (sampling period and number of samples) but halving the signal frequency to f=25MHz, see second figure

Can anyone explain to me why this happens? From Nyquist theorem I should be able to see the 50MHz component if sampling at 100MHz

Thanks for helping


texts: Findin the source direction using cross-correlation function

texts:
I am new to DSP and I'm working on the following problem.

We have a sound source in the far ﬁeld and a microphone array consisting of two microphones with distance d. We want to determine the source direction by using the cross-correlation of the microphone signals x and y.

I was told that in order to find the angle, we first need to calculte the time delay (in samples) between the microphone signals x and y. I know we can do this by finding the maximum of cross-correlation function R_XY. However, I don't understand why the global maximum of R_XY happenes at this time delay, namely k0.

Let's say we find this time delay. Then how do we calculate the angle? (assuming we know the distance between these two microphones, sound-propagation velocity and the sampling frequency.


texts: DFT of the same signal with different values of N

texts:
Let $x[n]$ be a discrete signal of 2 samples. We know that its DFT with N=4 is $X[k]=[0, 1+j, 2, 1-j]$. Without calculating $x[n]$, how can we know the DFT with N=2?

I have tried to use the relation between DFT and DTFT. Specifically, that the DFT is the result of sampling the DTFT in the frequency domain every $\frac{2\pi k}{N}$ samples. However, I don't know how to continue. I would really appreciate some help.


texts: In which area the image dithering technique is used?

texts: I know its an old technique but can someone tell me in which areas the image dithering technique is used and why?


texts:
Glottal Closure Instants (GCIs) are referred to the instances of signiﬁcant excitation of the vocal tract.

I have seen in some papers that I can measure instant of glottal closure in a voiced segment directly from speech waveforms.

Can someone provide me information how detect GCIs ?

How Glottal Closure Instants can be estimated ?


texts:
I have AM demodulated voice samples that come from a 14bit ADC that I will filter, process, and then send to a CODEC.  I can choose a sample rate of 96kSPS, 64kSPS, 48kSPS, or 32kSPS at the CODEC.  Taking into account most of the voice spectrum is limited to less than 15kHz, would I gain anything by sampling at 96kSPS vs 32kSPS?  I will use either an IIR or an FIR filter between the ADC and CODEC.

For a lower sample rate, the digital filtering would require more taps to acheive the same performance; however, there would be less samples to filter so that trade off seems like a wash.


texts:
I'm studying the transfer function of an aircraft that relates the pitch angle, $\theta$, with the horizontal stabilizer angle, $\delta$.
$$ H(s) = \dfrac{\theta (s)}{\delta (s)}$$

I have computed the frequency response by using Matlab, and I plotted the Bode Diagram: amplitude plot and phase plot.

I understand what is the amplitude curve indicating and I'm able to comment the result obtained. But when it comes to studying the phase curve (second figure above) I don't know what is the meaning and how it is related with the other curve.

So, my question is: what can I say from the phase plot by looking at the Bode Diagram of the studied transfer function, $H(s)$?


texts: Applying estimated transfer function to a signal

texts:
This is the simple code to find transfer function between sigout and sigin signals and then are the filter coefficients estimated by invfreqz function.

%% load simlated data
load CFDvsABAQUS_TC.mat

%% simulated data preprocessing (grid orientation B)
time = TCb(:,1);
sigin = TCb(:,2);           % inlet temperature (TC signal by Abaqus)
sigout = TCb(:,3);          % outlet temperature (by cfd FENIX)
Fs = 1/mean(diff(time));    % sampling frequency

%% Transfer Function Estimate outlet->inlet
[toi,foi]=tfestimate(sigout,sigin,2^13,[],[],Fs);

%% filter design by known transfer function
w = linspace(0,pi,length(fio)); % angular frequency f=[0,fs/2] -> w=[0,pi]
m = 500; % b coeffs
n = 500; % a coeffs
[boi,aoi] = invfreqz(toi,w,n,m);

%% final Digital filtering
esigin = filter(boi,aoi,sigout);

The input data are available at this dropbox link.

The final filter function produce completely wrong results (the result should be similar to original sigin signal). I have no idea what is wrong.

Any help will be very useful for me.


texts: How to identify and differentiate frequency and time in EEG data using python?

texts:
I have an EEG data of 200 Hz and sampled at 4097. I have few doubts and questions related to this data:

What does it mean by sampling in EEG data?
I am getting the sample rate and the frequency when I plot the data as shown below. Shouldn't it be the time and frequency?

I don't understand what I'm missing here. Can anyone explain how am I supposed to identify the time and frequency of EEG data?

EEG data file : http://filesave.me/file/123174/eeg-txt.html

Thank you in advance.

Code for the above plot:

x = np.array([])
input = open('eeg.txt', 'r')
for file in input:
    file = file.replace(',','')
    x = np.append(x, float(file))
input.close()
plt.subplot(1, 1, 1)
plt.plot(x)
plt.ylabel('EEG signal')
plt.show()

FFT calculation:

import numpy as np
x = np.array([])
input = open('eeg.txt', 'r')
for file in input:
    file = file.replace(',','')
    x = np.append(x, float(file))
input.close()
plt.subplot(1, 1, 1)
ff = fft(x)
plt.plot(ff)
plt.ylabel('EEG signal')
plt.show()

Wavelet decomposition:

coeffs = wavedec(x, 'db4', level=6)
cA2, cD1, cD2,cD3,cD4,cD5,cD6 = coeffs

plt.subplot(7, 1, 1)
plt.plot(x)
plt.ylabel('Noisy Signal')
plt.subplot(7, 1, 2)
plt.plot(cD6)
plt.ylabel('noisy')
plt.subplot(7,1,3)
plt.plot(cD5)
plt.ylabel("gamma")
plt.subplot(7,1,4)
plt.plot(cD4)
plt.ylabel("beta")
plt.subplot(7,1,5)
plt.plot(cD3)
plt.ylabel("alpha")
plt.subplot(7,1,6)
plt.plot(cD2)
plt.ylabel("theta")
plt.subplot(7,1,7)
plt.plot(cD1)
plt.ylabel("delta")
plt.draw()
plt.show()

Results of wavelet:


texts: quantitative definition of "FFTable" data point as a function of the signal to noise level

texts:
Regarding hotpaw2's answer to the "Number of FFT points required for a specific frequency resolution for an oversampled signal" question:

Given that the result of an FFT can be interpolated (possibly very accurately using Sinc interpolation), the number of FFT points required to estimate the frequency depends on the signal-to-noise ratio of the data containing your signal, and the type of resolution you require (peak separation, or peak estimation).
In the extreme case of zero noise or other interference, only 3 or 4 non-aliased points are required to exactly reconstruct a pure sinusoid and thus estimate its frequency.

I'm just wondering if you have a reference I can look up for the details about the quantitative relationship between the limit of signal to noise level and allowed data points for the FFT. That is, if I have a data set of 10 which has a signal to noise (S/N) level of 100, is it good enough to reconstruct 3 pure sinusoid and its frequency? and how about the case of S/N level of 10?


texts: Relation between power of the baseband signal and the low pass equivalent

texts:
Let $x(t)$ be a real bandpass signal and $x_l(t) = (x(t) + j\hat{x}(t))e^{-j2\pi f_0t}$ be the lowpass equivalent. Is there any relation between $P_x$ and $P_{x_l}$ where $P_x$ refers to the average power of $x(t)$ and $P_{x_l}$ refers to the average power of $x_l(t)$? I tried to use the definition $$P_x=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}|x_l(t)|^2dt = \lim_{T\rightarrow\infty}\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}(x^2(t) + \hat{x}^2(t))dt \tag{1}$$
And I got stuck here. Maybe we should add some assumptions on $x(t)$ in order to be able to simplify $(1)$.


texts:
I have the generic problem of analysing a discretly sampled signal, lets
say at 100Hz, and would like to compute features (ML application) of this signal at different scales:

An instantiation of this problem is computing the FFT on real-signal with
N samples:

Original signal: X[1],...,X[N]

Suppose I have already computed the FFT of the following
subarrays:

X[N-k+1],...,X[N]

X[N-2*k+1],...,X[N-k]

X[N-3*k+1],...,X[N-2*k]

X[N-4*k+1],...,X[N-3*k]

which is each of size k. Then I would like to compute the FFT
of the subarray

X[N-4*k+1],...,X[N]

(which is of size 4*k), WITHOUT reconstructing the original array X[.] and
running an ordinary FFT on the array of size 4k. Roughly I would like to do O(N) primitive operations in the best case if that is at all possible.

Is there some recursive property I could use?

PS: Note that if the goal of computing the FFT is substituted with e.g. finding the mean or the variance there exist easy recursive properties that allow merging of the results in constant time.


texts: Fourier transform of $\sum_{n=-\infty }^\infty e^{-|t-2n|}$

texts:
I need to find Fourier transform of

$$\sum_{n=-\infty }^\infty  e^{-\lvert t-2n\rvert }$$

I found the solution online. But I think the choice of limits in equation (1) of the solution is wrong. It should be probably from $-\infty$ to $2n$ for the first integral and from $2n$ to $\infty$ for the second integral. But I am confused about what the limits would be as the both integration and summation is involved. and I am not sure whose limits to change (probably limits of integration must be changed only and that of summation must be kept as it is)

Given Solution:

Please tell me what the limits will be and the explanation. Also if possible please solve the problem.

Original answer:

$$ \frac{(1-e^{-2(1+jw)})}{(1-e^{-2})(1+jw)}-\frac{(e^{-2})(1-e^{-2(1+jw)})}{(1-e^{-2})(1-jw)}$$


texts: Why does this transfer function has a second zero

texts:
I'm learning about $\mathcal Z$-transforms in DSP and I have a transfer function of the following form:

$$H(z)=\frac{2-3z^{-1}}{1-1.6z^{-1}+0.8z^{-2}}$$

When I calculate zeros and poles of this system by hand,
I get these poles from the equation

\begin{align}
\frac{z^2-1.6z+0.8}{z^2} &= 0\\
p_1&=0.8-0.4i\\
p_2&=0.8+0.4i
\end{align}

And a single zero from the equation

\begin{align}
\frac{2z-3}{z}&=0\\
z_1&=1.5+0i
\end{align}

However, when I use Wolfram Alpha to compute the poles and zeros, it also lists a second zero positioned at origin. Or in MATLAB:

The question is, where does that second zero come from?


texts:
I'm trying to implement color-to-grayscale method from this paper. And they use Verlet's integration as:

$$L^*(t+\Delta t)=\frac{F(t)}{m}\Delta t^2+2L^*(t)-L^*(t-\Delta t),$$

for computing new value of luminance. But they not say how to compute new value of luminance for $t=0$. According to this Wikipedia link I think it could be as:

$$L^*(\Delta t)=\frac{F(0)}{2m}\Delta t^2+L^*(0)$$

Is it correct?


texts:
It is known that the ACF of a 2D digital image can be obtained by computing the inverse Fourier transform (FT) of the product of the FT of the image and its complex conjugate. That is:

My problem is that I cannot figure out the expression for the inverse of G(p,q). I’d appreciate any idea, help, or suggestion. Thank you.


texts: Volume increase of a magnitude spectrogram

texts:
Let's say I have an audio waveform and I perform an STFT on the waveform and then take the absolute value to get the magnitude spectrogram.

Is there a transformation that I can make to that spectrogram that would produce the same result as taking increasing the volume of the original waveform and recomputing a magnitude spectrogram?

Is it as simple as multiplying all the bins by a constant?


texts:
I'm trying to do autocorrelation with numpy in python but i'm getting strange results :

For a simple 440Hz sine wave sample :

This is the result of autocorrelation :

Is this normal ? Maybe I'm not understanding everything

UPDATE :

Here is my code :

def autocorrelation(signal):
  res = np.correlate(signal, signal, mode='same')
  res = res[len(res)//2:]
  plt.plot(signal)
  plt.show()
  plt.plot(res)

  plt.show()


texts:
I need help solving the following blurring function question.Assume an image $f(x,y)$ is moving in front of a camera so that $_0()$ and $_0()$ are the time-varying components of motion in the x- and y- directions respectively. The camera’s shutter opens at  = 0 and closes at  = . Assume that the shutter opening and closing operations are instantaneous. The blurred image (, ) captured by the camera is calculated as:

$$
 g(x,y) = \int_{0}^{T} f( x - x_0(t), y - y_0(t)) dt
$$

Calculate the blurring transfer function $H(u,v) = G(u,v) / F(u, v)$ in  frequency domain with respect to $_0()$ and $_0()$, where $(,)$ and $(,)$ are the Fourier transforms of $(,)$ and $(,)$ respectively and then calculate the blurring transfer function for the case where the image $(, )$ moves in only x- direction with a constant speed $\frac{a}{T}$.

Attempt: I really have no idea how to set up this blurring function.It is clear that the speed, $v$, is

$$
 v = x_0(t)i + y_0(t)j
$$

and also that

$$
 \mathcal F \{ f(x - x_0, y - y_0) \} = \text{exp}( (-2\pi i (x_0u + y_0v) ) )F(u,v)
$$

$$
G(u,v) = \mathcal F \{ g(x,y) \} = \mathcal F \left\{\int_{0}^{T} f( x - x_0(t), y - y_0(t))\right\}
$$

but I don't know how to proceed from there.Any help appreciated


texts: Correct color values on image taken with no infrared filter

texts:
I have a Raspberry Pi NoIR Camera. Which I want to take pictures with at daytime in a well lit room.

Unfortunately the colors on the pictures do not look anything like you would expect from a normal camera. I guess the reason is:

  (NoIR = No Infrared.) This means that pictures you take by daylight
  will look decidedly curious, but it gives you the ability to see in
  the dark with infrared lighting.

This post states that you can only overcome this issue during daylight by using an additional infrared filter or by using a "normal" camera.

So the pictures seem to lack any green values, plus some of the black values appear red/purple.

I could apply filters to the images programmatically or using gimp and I was wondering how would one tackle such a task? Could this be solved with some simple color correction algorithms? Is it even possible to revert this effect without using a additional "hardware"-infra red filter?

Thanks in advance


texts: application of wavelet to vibration signal

texts:
I ask you to help me to interpret this figure corresponding to a diagnosis of gear failure using CWT (Haar);
I have already worked on the analysis of the vibratory signals using the FFT; But since the signals from the bearings and the gears are non-stationary it is considered necessary to apply WT. So I started with simulated signals of gear failures to evaluate the sensitivity of this method compared to that of Fourrier


texts: given a mixture of two signals, how to separate it to its base components?

texts:
I have a signal that contains two components. How can I separate it two its two components in matlab?
Please note that the location and shape of components are different in each experiment and we do a lot of experiments. It should automatically detect the two components based on the frequency contents. I have added this comment because some people suggested
c1 = x(1:80);c2 = x(81:end);

I do not mean such separation.


texts: Calculating SPL from pressure signal - Amplitude vs Power method

texts:
I have a pressure signal from a Fluent FFowcs-Williams Hawkings acoustics analysis.

I converted this pressure signal into the frequency domain in order to get SPL values, using Matlab. I used the amplitude of the FFT:

A = fft(s10);        %Amplitude
B = abs(A);          %Magnitude of amplitude (complex numbers to real)
C = B/L;             %Normalize amplitude by dividing by number of samples
D = C(1:fix(L/2)+1); %Select one half of 2-sided spectrum
E = 2*D;             %Compensate for loss of energy by taking only 1 half of spectrum
F = E/sqrt(2);       %Rms of amplitude

SPL = 20*log10(F/Pref);

As far as I can see, Fluent uses the Power Spectral Density of the FFT to calculate SPL, instead of using the amplitude. The results of this calculation corresponds to what Fluent does in its FFT:

A = fft(s10);                 %Amplitude
G = A(1:fix(L/2)+1);          %Select 1-sided spectrum
PSD = (2*(abs(G)).^2)/(Fs*L); %Compensate for loss of energy by multiplying by 2; scale by
                              %dividing by Fs and L; square
SPL2 = 10*log10(PSD/(2e-5)^2);

When plotting these two different sets of results, there is a constant difference between the SPL values of the two signals. In this case it's a 12.0377 dB shift.

My question is, which of these two methods is the correct one to use in order to calculate SPL, and why? Or have I made a mistake in my calculations?

Help would be much appreciated!!


texts: Discontinuity at the edge of Chebychev window

texts:
I am using Chebychev window for its narrow main lobe. The problem of chebychev window is that it has discontinuities at the edge, and it seems that Taylor window solves this issue.

More detail:
http://de.mathworks.com/help/signal/ref/taylorwin.html
http://en.wikipedia.org/wiki/Window_function#Dolph.E2.80.93Chebyshev_window

I've searched around but I can't find any information on how to implement a Taylor window. Any information on taylor window or suggestions on fixing this issue of edge discontinuities would be very appreciated.


texts:
Forgive me if it’s too basic, I finish engineering a while ago.
Given any time series, not periodic, I would like to find any repetitive pattern that is distinct (by some given measurement) and is unknown.

For example given this :
 I would like to find number 1 and 2( get their index on time).

Is FFT going to help anyway?
Without knowing or hard programming an algorithm, is it possible or do I need a serious machine learning?

(Are there practical books about finding patterns with code and not too much math?)
Thanks.


texts: Frequency Domain Filtering with big kernel size

texts:
I use FFT to do filtering in frequency domain, but when I use big kernel I got shift on border, I think it's due to nature of cyclical convolution and I need to do more zero padding.

But how to calculate required padding size? should it be power of 2 in matlab?

Here is my matlab code:

tic
I= rgb2gray(imread('im.jpg'));
h = fspecial('gaussian',555);
F = fft2(double(I));
H = fft2(double(h), size(I,1), size(I,2));
F_fH = H.*F;
ffi = ifft2(F_fH);
toc

figure, imshow(ffi,[])

UPDATE:

Here is solution:

%The zero-padding is just a way to make cyclical convolutions
%(which is what FFT-based convolutions are) act like linear convolutions.
%the convolution theorem yields the desired linear convolution result only
%if x(n) and h(n) are padded with zeros prior to the DFT such that their
%respective lengths are Nx+Nh-1, essentially zeroing out all circular artifacts.

tic
%Frequency Domain Filtering
I= rgb2gray(imread('im.jpg'));

ker_sz= 55;
h = fspecial('gaussian',ker_sz);
F = fft2(double(I), size(I,1)+ker_sz-1, size(I,2)+ker_sz-1);
H = fft2(double(h), size(I,1)+ker_sz-1, size(I,2)+ker_sz-1);
F_fH = F.*H;
ffi = ifft2(F_fH);
ffi= ffi((ker_sz-1)/2+1:(ker_sz-1)/2+size(I,1), (ker_sz-1)/2+1:(ker_sz-1)/2+size(I,2));
toc
%Display results
figure, imshow(ffi,[])


texts: Expected value of coupled signals is zero

texts: Suppose there are two vector signals, $u$ and $v$, which are functions of time. Given, $E[uv^T]=0$, what is the physical significance of this statement. For clarity, assume $u$ and $v$ are noises in a system like in kalman filter, control noise and measurement noise. E is the expected value and assume that the vector dimensions are entirely compatible.


texts: What is difference between terms $X(j \omega) ,X(\ e^{j \omega })$ and $ X(\omega)$?

texts:
While studying frequency transforms ,I get confused with the terms like $X(j \omega) ,X(\ e^{j \omega })$ and $ X(\omega)$
 ,where $ \omega = 2 \pi f $.
 So what is the difference between them ?


texts: About the proof of an equality related to the DFT [sampling the DTFT to obtain the DFT]

texts:
This wiki page about the DTFT says that the DFT can be obtained from the DTFT by sampling the latter in one cycle at $N$ points:

  When the DTFT is continuous, a common practice is to compute an arbitrary number of samples (${{N}}$) of one cycle of the periodic function (${X_{1/T}}$):\begin{align}
\underbrace{X_{1/T}\left(\frac{k}{NT}\right)}_{X_k} &= \sum_{n=-\infty}^\infty x[n]\cdot e^{-i 2\pi \frac{kn}{N}} \quad \quad k = 0, \dots, N-1 \\
&= \underbrace{\sum_{N} x_{_N}[n]\cdot e^{-i 2\pi \frac{kn}{N}},}_{DFT}\quad \scriptstyle{\text{(sum over any }n\text{-sequence of length }N)}
\end{align}
  where $x_{_N}$ is a periodic summation:
  $$x_{_N}[n]\ \triangleq\ \sum_{m=-\infty}^{\infty} x[n-mN]$$

So we essentially have the following equality
$$\sum_{n=-\infty}^\infty x[n]\cdot e^{-i 2\pi \frac{kn}{N}}= \sum_{n=0}^{N-1} \sum_{m=-\infty}^{\infty} x[n-mN] e^{-i 2\pi \frac{kn}{N}}. $$

I didn't manage to prove it, is this a known identity?


texts: Filter size vs. FFT size and Overlap add

texts:
First of all I would like to say that :

I'm new to DSP and in a "learning curve" hence the "basicness" of my questions.
I've read as many posts/articles as possible on the internet but still need "guidance".
I'm french hence the bad formulation sometimes.

My mission/data context :
Performing data analysis on real data sets (oil-rig data).
Sometimes data collected on the rig present low frequency component due to rig movements.
Data are collected every 5 seconds so the sampling rate is 0.2 Hz. The size of the data to be analysed will be 1024/2048/4096 points let's call it N is this post.
My mission is to perform a filtering on these data to eliminate the low frequency component.
This is not a real time application. The range of data to be analysed will be chosen and then the analysis performed.
Chosen approach :
Process these data with a FFT, filtering the FFT and perform an IFFT to obtain and clean(er) set of data. Basically the lowest frequency found in the FFT is to be eliminated.
Many users suggested me to use the overlap-add method to obtain the desired result.
First I was thinking of just zeroing low frequency related bins in my FFT but my readings on the subject suggested to rather use a high pass filtering (with spectral inverted windowed-sinc).
Designing this filter and applying it is what I post here for because it turned out to be way more difficult than expected.
Questions :
Do I have to use the overlap-add method considering that my application is not a real time one ?
For what I have understood I need to design a high pass filter kernel with a spectral inverted sinc function. Then I have to perform a FFT on my filter kernel and multiply points by points my filter FFT(1) coefficients with my data FFT(2) coefficient before performing a IFFT on the resulting FFT. Am I correct on this one ?

According to my readings I can obtain a spectral inverted sinc function by inverting (-sin(x)/x) it and adding 1 to the center of symmetry. I am right, do you guys confirm ? Do I add 1 on the one center point only ?

What is the right size to choose for my filter kernel ? Is it the same size N as my data ? As I have to multiply my FFT(1) with my FFT(2) I would think that I have to choose the same size but I'm probably wrong on this one and the overlap add method may be the answer.

The kernel filter designing is tricky, I don't know where to place the central lobe of my spectral inverted sinc, in 0, on a specific offset in the time domain ?


texts: Computing Hartley transform of zero-padded signal from transform of original signal?

texts: If I have the discrete Hartley transform of a signal of length 2n, how can I use it to efficiently obtain the discrete Hartley transform of the same signal padded to length 2n+m for some m > 0?


texts: Degrees of freedom in coding (for a communication channel)

texts:
I was reading on my lecture slides about how to achieve diversity (i.e. transmit more copies of the same message in order to avoid its corruption; this technique may be realized by repeating the message in different instants, by sending it at different frequencies, and by using multiple antesas) in a communication system.

This is what it is written about this topic:

  Consider now a MIMO channel with two transmit and two receive
  antennas:

  One can also design codes specifically for the transmit diversity
  system (space‐time codes). The simplest of these, and yet one of the
  most elegant, is the so‐called Alamouti scheme. This is the transmit
  diversity scheme proposed in several third‐generation cellular
  standards and is designed for two antennas (but can be extended to
  more antennas to some extent):

  With BPSK symbols, it can be shown that the repetition scheme requires
  a factor of 2.5 (4 dB) more power than the Alamouti scheme. The
  repetition scheme suffers from an inefficient utilization of the
  available degrees of freedom in the channel: over the two symbol
  times, bits are packed into only one dimension of the received signal
  space, namely along the direction .
  In contrast, the Alamouti scheme spreads the information
  onto two dimensions – along the orthogonal directions

  The degrees of freedom (DoF) of a channel can be defined as the
  dimension of the received signal space. In a channel with two transmit
  and a single receive antenna, this is equal to one for every symbol
  time. The repetition scheme utilizes only half a degree of freedom per
  symbol time, while the Alamouti scheme utilizes all of it.

Now, I need a basic and simple explanation of what is happening. Precisely, what does it mean in practice "degrees of freedom" of a channel?

Consider for instance the Alamouti scheme, in which the message u is described through a 2x2 matrix (with columns (u1, u2) and (-u2*, u1*). Degrees of freedom is the rank of this matrix, if I have understood it correctly. But what does it mean in practice? Which properties should the message have to achieve this?


texts: Scipy Welch's gives different first element

texts:
I've re-implemented Welch's method and want to compare it to scipy.signal.welch. However, the first two and last elements of the resulting array are different.

My Welchs function (not optimized and without taper, as this is not the issue here)

import numpy as np
import scipy

def welchs(signal, nperseg=512, noverlap=0):
    taper = np.ones(nperseg) # assume boxcar taper for simplicity
    scale = 1.0 / (1 * (taper*taper ).sum()/2)
    siglen = signal.shape[-1]
    step = nperseg - noverlap
    fft_segments = []
    for i in range(0, siglen, step):
        if i+nperseg>siglen:
            print(i+nperseg, '>', siglen)
            break
        segment = signal[i:i+nperseg]
        w = np.fft.rfft(segment)
        fft_segments.append(np.abs(w)**2)
    fft_segments = np.array(fft_segments)
    fft_mean = np.mean(fft_segments, 0)
    return fft_mean * scale

testing:

signal = np.sin(np.linspace(0,512, 1024))

w1 = welchs(signal, nperseg=512, noverlap=0)
_,w2 = scipy.signal.welch(signal, nperseg=512, noverlap=0, detrend=False, window='boxcar')
np.testing.assert_almost_equal(w1,w2) # -> 0.7% mismatch, 2 elements

All elements of w1 and w2 are the same except the first and the last.

Can anyone explain to me what my implementation might be missing?


texts: Design a linear-phase FIR filter approximating the magnitude of a given IIR filter

texts:
I have a biquad IIR filter from which I want to get a linear-phase FIR. I saw that related question but the OP does not care much about phase.

I have troubles getting the exact same magnitude curve when extracting an FIR from my IIR.

So far, my method has been:

Apply my IIR filter on a 1024 samples buffer which contains nothing but a Dirac as first sample.
Shift the impulse response obtained in step 1 and make it symmetric (to get phase linearity)

If I stop after step 1, I get the exact same magnitude curve for my IIR and my FIR but I also get the exact same phase curve (which is non linear, and hence not interesting).
If I stop after step 2, I get an linear phase but not the exact same magnitude.

On the following picture, IIR's and FIR's (after step 2) transfer functions are (respectively) in red and blue:

What have I been doing wrong?


texts: Find horizontally and vertically aligned sets of points

texts:
I have the image below.  This image contains sets of points that are aligned horizontally and vertically.

How can i find these sets of points? based 1) on their alignment and  2) on the ** **inter-line and Inter-word distances.

Do I have to use these 2 arguments, alignment and distance? or is one of them sufficient?

Is it possible to use the Hough transform? or another algorithm?


texts:
I would really appreciate if I can get answers to following Fixed point queries(part-1):

Converting a real number x to fixed point fraction X(m,n)

[1] says:
m = ceil(log2(max[abs(xmin,xmax)+1])) { assuming signed x>1) or
m = floor(log2(max[abs(xmin,xmax)+2])) { assuming signed x<1}            (1.1)
and xmin,xmax define the known range of x. if WL desired is known = m+ n , then, n= WL - m.

Q1) Can this be used blindly? What if we don't know the expected range of x?

[2] defines:  if floating x had "d" fractional digits,
compute ceil(3.3 * d).                   (1.2)
Q2)  Is Fraction Length n = ceil(3.3 * d)  or n+1  = ceil(3.3 * d)
Q3) What about floor(3.3*f) or any other scheme?
Q4) Where would I use equation 1.1 over 1.2 and vice versa?
Suppose we have computed Fraction  length n by any of the above methods.
Now, The integer representation and the quantized value is computed as:
X= round/floor/fix(x*(2^n)); xq(Quantized value) =  X/(2^n);
(Q5) Does choosing round/floor/fix here determine the type of quantizer? If yes,
   which option gives a rounding quantizer and which one does a truncation quantizer?

Reference:

1.http://www.superkits.net/whitepapers/Fixed%20Point%20Representation%20&%20Fractional%20Math.pdf( equation 9 and Equation 12)
 2.  DSP, a practical approach:Efeachor (page 809-810,second edition)

Thank You so much
sedy


texts:
I have a sinusoidal signal embedded in noise, and I would like to estimate it's amplitude and frequency, so I can subtract it out.

I realize that an FFT can do the job for me, and that's what I've been doing - I took an FFT of the input signal, and looked at the frequency bin with the highest amplitude, used that as an estimate for the frequency of the wave, and then subtract it out.

But the problem is that the frequency is accurate to only half a bin, and the bin size can be fairly large (since my sample size is typically ~ 100 points), and hence the error in wavelength is big. This leaves fairly large residuals once I subtract.

So as step two, I used the frequency from the FFT as an estimate, and I used the Levenberg-Marquardt algorithm to fit for the frequency, phase and amplitude. This solution worked wonderfully for toy data that I generated, but it fails quite frequently on real data, because of the noise. The amplitude of the noise can be as large (at most) as the amplitude of the signal.

I figured that I can estimate the wavelength using an autocorrelation, but I'm not sure how to interpret the results in the case of more than one frequency component.

Question: Are there any methods of estimating the wavelength of the wave that aren't FFT? I figure since fitting for the frequency is failing with noise, I could try fitting for the wavelength, but any other suggestions are welcome!

Note: I haven't really tried padding with zeros while performing the FFT because the noise seems to go up in the fourier plane if I do, but I'll try it and report the success here.


texts: How to Generate the Following Image Effect (Modulation)?

texts:
I was given the original image and I need to produce Image 1. Does anybody know how I would go about this? I am on MATLAB, but I am more interested in the theory behind this.

Original Image

**Processed Image **

NOTE How would you achieve this through convolution of the ffts of the image and a sine wave along the vertical axis?


texts: Can two different impulse responses give the same frequency response?

texts: I have two different impulse response each with different length, is it possible that they have the same frequency response?


texts: Correct Fish Eye camera Radial Distortion and blurry effect at the image edges

texts:
There are a lot of questions regarding this topic ,but I am trying to get a more clear picture from these questions.

I am trying to calibrate a fish eye camera and I am using OpenCV omnidir class functions to find the camera intrinsics.

I am getting fair results. The problem is ,

At the image edges the objects gets stretched and also I am losing some information at the edges.

Here is my input image:

Here is my output image:

As you can see, I am losing some information at the edges(left and right) and also the images start stretching at the edges.

My questions are as follows:

How I can I include more FOV in the corrected image at the edges, where the information is lost?
How can I reduce the blur effect at the edges?
During calibration, should I cover the entire FOV of the camera so that the corners are present at the edges also?
What is the correct way of showing the patterns while calibration?
Are there any online tool boxes which provides fish eye calibration.

Here is my code snippet for calibration and testing

//Calibration
Mat K, xi, D, idx;
int flags=0|omnidir::CALIB_FIX_SKEW | omnidir::CALIB_FIX_K1 |
 omnidir::CALIB_FIX_K2;

TermCriteria critia(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 200,
0.0001);

vector<cv::Mat> rvecs, tvecs;

double rms = cv::omnidir::calibrate(obj_points, image_points, image_size, K,
xi, D, rvecs, tvecs, flags, critia, idx);

//Testing
Mat R = Mat::eye(3, 3, CV_32F);
Mat Mapx, Mapy;
Mat New_camera_mat(3,3,CV_32F);

//New_camera_mat tries to get entire FOV,but it is losing some information
  at edges

New_camera_mat.at<float>(0, 0) = 100; New_camera_mat.at<float>(0, 1) = 0;
New_camera_mat.at<float>(0, 2) = 1280/2;

New_camera_mat.at<float>(1, 0) = 0; New_camera_mat.at<float>(1, 1) = 100;
New_camera_mat.at<float>(1, 2) = 720/2 ;

New_camera_mat.at<float>(2, 0) = 0; New_camera_mat.at<float>(2, 1) = 0;
New_camera_mat.at<float>(2, 2) = 1;

cv::omnidir::initUndistortRectifyMap(K, D, xi_Right, R, New_camera_mat,
image_size, CV_32F, Mapx, Mapy, cv::omnidir::RECTIFY_PERSPECTIVE);

remap(distorted_frame, undistorted_out_frame, Mapx, Mapy, INTER_CUBIC);


texts: Appropriate digital filter for GPS tracks?

texts:
I have an open-source software project whose purpose is to analyze a GPS track, or a similar track made by an application such as google maps, and estimate the physical exertion required to hike or run that route. Traditionally, people have gotten a gut feeling for this sort of thing by specifying the horizontal distance and the vertical gain (i.e., the amount of climbing, ignoring any descents). My software tries to do better than that, but in any case, all such measures are extraordinarily sensitive to small, spurious fluctuations in the elevation data. For instance, I could carry a GPS receiver with me on a totally flat run along a straight road in Kansas, and if the elevation is constantly fluctuation up and down by a few meters, due to GPS errors, it could show up as if I have a huge amount of vertical gain. Similar things happen with map-tracing methods, because the digital elevation models (DEMs) have fairly large errors. GPS elevations can be extremely accurate when the position of the satellites in the sky is favorable, but it can be incredibly far off in unfavorable cases (I've seen errors of thousands of meters).

Can anyone recommend a suitable method for filtering, so that I'm not reinventing the wheel, and doing it badly?

The crude method I've been doing so far is the following. I do an initial iteration of the data analysis to find the integrated horizontal distance $h$ for each point on the track. That gives me the coordinates $x(h)$, $y(h)$, and $z(h)$ as functions of distance. Then I convolve the elevation $z(h)$ with a rectangular window 500 meters wide horizontally. This gives me what appear to be accurate and reproducible estimates of the total gain. However, I also find that my initial estimate of $h$ is quite crude, because the track has a sort of fractal structure, some of which may be real and some of which is probably measurement errors. So I probably want to do some kind of low-pass filtering on $x$ and $y$ as well, maybe cutting out anything with a wavelength less than about 100 m, and then do a second iteration on the integration of $h$.

Googling turns up a lot of material on Kalman filters. Am I correct in understanding that a Kalman filter is not really the right tool for this job, since it's meant for an object like a missile or a helicopter, which has a lot of inertia? Also, my data is a track in Euclidean space, not a time series.

I'm looking for a method that is fairly robust, and can be easily implemented and played with using open-source software on linux. My code is written in ruby, but I would be OK with shelling out to something like scipy, as long as it has a fast startup time. (E.g., I don't want to use Julia because of the slow startup time.) I would like to implement this using a library that is well tested, has a large user base, and is likely to be around and well maintained for a long time.


texts: FFT Convolution: missing sample somewhere

texts:
For a piece of software I've been working on lately I need to implement FFT convolution for performance reasons.
Doing so has been mildly succesfull but there is still a error somewhere causing lots of small spikes to show up on the waveform. Also when looking at the spectrum of the signal, lots of distortion can be seen which is filtered out mostly at 15 KHz.
What I do in the filter is the following with overlap - add:

Create the filter kernel of 101 points and FFT it
Create a buffer holding 200 samples of audio
When there are 200 samples ready, FFT the signal and multiply with the filter kernel FFT.
IFFT the filtered signal again
Add last 100 samples from the previous run to the first 100 samples of the current output
Save the last 100 samples for the next run.

The small spikes and distortion when looking at the spectrum seems to come from a missing sample somewhere, but I can't track down exactly where it happens.
The following code is what I use to do the filtering:
function TFFTAnalyzer.doAnalyzer(Sample: Single) : Single;
var
  I: Integer;

  Temp : Single;
begin
  Result := X[BufferPos];
  X[BufferPos] := Sample;

  BufferPos := BufferPos + 1;
  if BufferPos > 200 then
    BufferPos := 1;

  if BufferPos = 200 then
  begin
    //Add 100 zeros
    for I := 0 to 512 do
    begin
      if I > 200 then
        X[I] := 0;
    end;

    //Calculate the FFT of the sample buffer
    FFT_transform(1024, X);

    I := 0;

    while I <= 1024 do
    begin
      Temp := X[I] * FilterX[I] - X[I + 1] * FilterX[I + 1];
      X[I + 1] := X[I] * FilterX[I + 1] + X[I + 1] * FilterX[I]; // imaginary part
      X[I] := Temp; // real part

      I := I + 2;
    end;

    //Inverse FFT the signal
    FFT_inverse(1024, X);

    //Overlap the first 100 samples
    for I := 0 to 100 do
    begin
      X[I] := (X[I] + Overlap[I]);// / 2;
    end;

    //Save last 100 samples for next overlap
    for I := 200 to 300 do
    begin
      Overlap[I - 200] := X[I];
    end;
  end;
end;

This is the picture showing the problem with the spikes and distortion.


texts: How to connect $\omega$ in time domain to physical frequency?

texts:
I'm a noobie to DSP, and recently I'm learning DSP for an urgent assignment. There's a problem that has been confused me for days. I hope anybody could shed some light on it, although it might seem to be silly.

The assignment asks me to filter a signal so that it's in $20-800$Hz. After I spent several days to learn DSP, I finally get some idea about frequency domain, DTFT, FIR, IIR, etc. But I still have no idea about how to proceed. Furthermore, I have a huge confusion about frequency: in the frequency domain, we always talk about the radian $\omega$ which is always less than $\pi$, how can I connect this $\omega$ to the physical frequency? Is it right to use $\omega/2\pi$?


texts:
The system:
 $$ T(x[n]) = ax[n] + bx[n-3] $$

For me it seems that the system is linear:

$$
\begin{align}
T(\alpha_1x_1[n]+\alpha_2x_2[n]) & = a(\alpha_1x_1[n]+\alpha_2x_2[n]) + b(\alpha_1x_1[n-3]+\alpha_2x_2[n-3]) \\
& = \alpha_1(ax_1[n] + bx_1[n-3]) + \alpha_2(ax_2[n] + bx_2[n-3]) \\
& = \alpha_1T(x_1[n])+\alpha_1T(x_2[n])
\end{align}$$

Thus it's linear, however in the presentation I got it says it's not linear. (without reasoning) Where I'm making the mistake (or maybe there is a mistake in the presentation)


texts: discrete fourier transform circular symmetry

texts:
I was reading digital signal processing. In topic Discrete Fourier Transform, circular symmetry it is said that circular advance is obtained by shifting x(n) in clockwise direction. i.e. to obtain x(n+k) Shift sequence x(n) in CW direction by k samples.

my question is would it same for x(-n+k)?

circular fold x(n) to obtain x(-n) and shift sequence x(-n) in CW direction by k samples. am i right?

*CW=Clock Wise


texts: Eigen Function of Linear Time Invariant (LTI) System

texts:
Whenever I try to find the derivation of eigen signal of a LTI system, they always start with assumption that the input signal is complex exponential. Is there any other way to prove that this is the case starting from the basic equation of LTI system
$$\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau = \lambda x(t)$$
where $\lambda$ is the eigenvalue of the eigen signal


texts: Fast Walsh-Hadamard Transform frequency spectrum?

texts:
Can the Fast Walsh-Hadamard Transform be used to find a Walsh-Hadamard frequency spectrum which some how relates to a Fourier frequency spectrum?

Is the Fast Walsh-Hadamard Transform more accurate than the Fast Fourier Transform when analyzing signals with very low bit resolution?

Below are the different Walsh-Hadamard harmonic signals stored in the matrix's columns ordered from lowest frequency to highest frequency across the matrix's rows with a lot of phase variation degeneracy in the mid-range frequencies. The matrix is symmetric.

octave:56> ifwht(eye(16)) .+ 1
ans =

   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2
   2   2   2   2   2   2   2   2   0   0   0   0   0   0   0   0
   2   2   2   2   0   0   0   0   0   0   0   0   2   2   2   2
   2   2   2   2   0   0   0   0   2   2   2   2   0   0   0   0
   2   2   0   0   0   0   2   2   2   2   0   0   0   0   2   2
   2   2   0   0   0   0   2   2   0   0   2   2   2   2   0   0
   2   2   0   0   2   2   0   0   0   0   2   2   0   0   2   2
   2   2   0   0   2   2   0   0   2   2   0   0   2   2   0   0
   2   0   0   2   2   0   0   2   2   0   0   2   2   0   0   2
   2   0   0   2   2   0   0   2   0   2   2   0   0   2   2   0
   2   0   0   2   0   2   2   0   0   2   2   0   2   0   0   2
   2   0   0   2   0   2   2   0   2   0   0   2   0   2   2   0
   2   0   2   0   0   2   0   2   2   0   2   0   0   2   0   2
   2   0   2   0   0   2   0   2   0   2   0   2   2   0   2   0
   2   0   2   0   2   0   2   0   0   2   0   2   0   2   0   2


texts: Approximate the magnitude response of an analog filter with a digital filter, starting from its gain expression

texts:
I have an analog filter with its frequency response curve in dB described by the following expression:
$$
N_{dB}=20log_{10}\omega t_1 \sqrt{\frac{1+(\omega t_2)^2}{1+(\omega t_1)^2}}
$$
This expression is derived from the series connection of two lowpass filters each associated with the following RC circuit:

where, for each circuit, the time constant $t_i=RC$ and of course $\omega = 2\pi f$, where $f$ is the frequency (this is actually the equalization curve for magnetic tape recording/playback, see Annex B, page 14 of this document).

I would like to obtain an approximation of this frequency response using a digital filter. I don't know if there is a method to exploit our knowledge of the analog frequency response or if I should design the filter myself from scratch.
The end goal is to obtain the impulse response and save it as a a .wav file (I know how to do this last part). I just took a basic DSP course at my Uni but we didn't work with analog filters so I am a little bit lost.


texts: I'm a newbie to Kalman filter. I have found the code online but I was wondering if there is any sample dataset available online to get hands-on with it (for example: CIFAR-10 for classification etc. ).


texts:
Consider the signals $e(t)$ and $h(t)$ with units $\left[\frac{V}{m}\right]$ and $\left[\frac{A}{m}\right]$ and their Fourier transforms $E(\omega)$ and $H(\omega)$ with units $\left[\frac{Vs}{m}\right]$ and $\left[\frac{As}{m}\right]$. The time avaverage P of the product of the two signals (~power spectral density), $e(t)h(t)$, in function of the Fourier transforms $E(\omega)$ and $H(\omega)$ can be written as
\begin{align}
P &= <P'(t)> = \lim\limits_{T\to \infty} \frac{1}{T} \int\limits_{-T}^{T} e(t)h(t) dt \\
&= \lim\limits_{T\to \infty} \frac{1}{T} \int\limits_{-T}^{T} \int\limits_{-\infty}^{\infty} E(\omega)e^{i\omega t} d\omega  \int\limits_{-\infty}^{\infty} H(\omega')e^{i\omega' t} d\omega' dt \\
&=  \int\limits_{-\infty}^{\infty}  \int\limits_{-\infty}^{\infty}E(\omega) H(\omega')\left[ \lim\limits_{T\to \infty} \frac{1}{T} \int\limits_{-T}^{T}e^{i(\omega+\omega')t} dt \right] d\omega d\omega'  \\
&=  2 \int\limits_{-\infty}^{\infty}  \int\limits_{-\infty}^{\infty}E(\omega) H(\omega')\delta(\omega'+\omega) d\omega d\omega' \\
&=  2 \int\limits_{-\infty}^{\infty}  E(\omega) H(-\omega)d\omega \\
\end{align}
The question: the unit of the time average $P$ should be $\left[\frac{VA}{m^2}\right]$ whereas the unit of the final result is $\left[\frac{VA}{m}s\right]$. How can this be and how to interpret it?
Example: if you know the Fourier transforms $E(\omega)$ and $H(\omega)$, and you want to calculate the average $P$, then you will use $P=2 \int\limits_{-\infty}^{\infty}  E(\omega) H(-\omega)d\omega$ but this will give another unit as expected?! So, then you have the wrong result?

What is unit of time average of product of two signals?

It's the product of the units of the signal. In my example is $V \cdot A = W$ which is power and in yours it's  $V/m \cdot A/m = W/m^2$ which is Intensity.


texts: Numerical higher order derivatives and time axis

texts:
I have a rather elementary question. Suppose we wish to study even-derivatives of an instrumental signal say second fourth and sixth derivatives and plot it as a function of time. With each successive differentiation we lose the data points. If the initial number of data points were 1000, the second derivative will have 9998 points, the fourth derivative will have 9996 points and so on. This can be seen with the diff command in Matlab.

How should we align these derivatives with the time axis which has 1000 points for plotting purposes? The typical way is to leave the first two points of the time axis for the second derivative and first four points for the fourth derivative. Is there a mathematical justification for this way of alignment since a true derivative is defined at a point rather than this numerical approximation? In some chemistry texts, authors suggest to average the values of the x-axis and plot the derivatives with respect to that.

Thanks.


texts: What are some techniques to recover an original image from three different compressed versions?

texts:
This question is theoretical. Let's say you have three different highly compressed images of a same source image, each compressed with a different algorithm and thus "lossy" in different ways.

What would be the possible techniques to recover the best possible image, using the three different compressed versions ?


texts: Simulation of a pressure field developed by a song

texts:
Greetings!

What I am doing:

In short, I want to write a program that is fed a song as a wav file, and spits out an animated contour plot indicating the relative pressure/loudness/quantity indicating amplitude at various points in a plane. Each speaker is modeled as a point source a set distance from each other hanging in an infinite void. My weapon of choice is Matlab.

How I am doing it:

Poorly. Song is read and split into channels, then cut into little slices. Each slice is FFT'd and the amplitudes and phases are split. For each point in the simulation space every frequency is taken, wavelength calculated and used in combination with the phase to determine the contribution of that frequency to the dimensionless loudness/pressure/sound intensity indicating quantity in that cell (the exact unit isn't that important - explained in next section). Do this for each slice of song, draw contour plots for each slice, and animate!

What my problem is:

The above seems kinda clumsy, and my preliminary results tell me I'm likely neglecting something huge or making some kind of fundamental error. Please note the idea here isn't to make a simulation for speaker setup or any kind of real physical use, I want to draw pretty pictures of what songs do with some basis in reality (specific application: psytrance). This is why the units are not that important, whether I use the normalised voltages in the wav file or convert it to dB or pressure... apart from looking a bit different in lin or log space the relative difference should be the same, at least for my pretty-picture-drawing-purposes. Once the program works as planned some number-fudging will occur to make it look better.

My Questions

1) Can someone suggest a better way of doing this? I feel like I'm doing this in a too complex manner, but I don't even know where to begin looking for a different method to accomplish my goals.

2) Am I right in doing phase corrections with the values given in the FFT? Gut says the error is here - Fourier coefficients are returned as complex numbers a+bi - for the use I have described here, how would I use this to yield "amplitude" and phase angle? It's a simple question... but humour me, self-kicking might follow your answer.

Edit: On second thought this may not be the most appropriate place to post this (at least for part 1 of the question), if someone can suggest a better place, please do.


texts: low-frequency amplitude with phase modulation of low Hz carrier

texts:
I am writing an FM/PM synthesizer, using Java. For the most part it works pretty well. What I am curious about is an attempt to use a low frequency tone as a carrier. By low frequency, I mean in the vicinity of 1 Hz. I use this a lot in my FM-synth coding as it warms up/thickens the output via the slow beating of overtones that occurs between sidebands.

The result I am getting sounds similar to the output of the Yamaha DX7S, and even closer to the output of the Native Instruments FM7 (I am behind an upgrade or two). My favorite sound, though, is the DX7S result. In trying to find out the nature of the differences, I recorded and examined output from all three with Audacity.

Below are three pairs of images from Audacity from three sources. In each case, the carrier is 1Hz and the modulator is 440Hz.
1) Yamaha DX7S
2) Native Instruments FM7
3) Java that I coded.

I tried to get the amplitudes and timbres to roughly match before making these recordings, though didn't go to a lot of trouble in that regard. (All of these have a sort of "wooo" lower format quality, somewhere between "boot" and "boat.") I was using a modulation factor of 200 in the Java version, where the size of the lookup table is 1024 units, containing one period in the range -1 to 1. The DX7S and FM7 both use a number in the range 0..100 for designating the modulation Index, but they don't seem to be defining them in the same way. The FM7 is set to "38" and the DX7S is set to "84".

1) Yamaha DX7S

There is noise from the unit itself and from the external mixing gear it is running through. But there is no prominent low Hz frequency component compared to the following diagrams.

2) FM7

3) My Java coded phase modulation

I like the sound of the DX7S effect the best. I'm noting that the FM7 also has a large low Hz component in the frequency analysis, and so am wondering if perhaps Yamaha was using some sort of high-pass filtering.

Below is the core Java code that I'm using for the phase modulation. It is being called once per audio frame (44100 fps), and makes use of a lookup table.

public float get(float modulationAmt)
{
    cursor += pitchIncr;

    // Handle sine LUT overshoot by wrapping around
    while (cursor >= tblSize ) cursor -= tblSize;

    // Phase modulation requires we do not directly modulate
    //    the cursor.
    temp = cursor + modulationAmt;

    while (temp < 0) temp += tblSize;
    while (temp >= tblSize) temp -= tblSize;

    return SineTable.get(temp);
}

Any thoughts as to what might be going on? The visual wave makes me wonder if there is DC bias of some sort happening, even though phase modulation is supposed to eliminate 0 Hz energy, if I read the Schottstadt article correctly (easy for me to misread, as it overshoots my math comprehension).

Is there an efficient way to eliminate the 1Hz component without removing the desired beating of sidebands? If Yamaha had accomplished this in the 1980's I'm thinking the algorithm must have been pretty minimal. For example, I have read that in some contexts, a simple "box" filter of 2 frames was used to help prevent aliasing. I haven't verified this though, and a box filter is low-pass, not high-pass. But maybe a box filter with a large number of frames can be inverted, or something like that, to make a "cpu cheap" high pass filter?

Or maybe that large low Hz peak can be ignored and isn't really relevant to the sound differences I'm hearing. But on general principles it seems to me having a large amount of "unheard" energy or output could potentially be asking for problems.


texts: How to remove a tone/interference from a signal, without using classical linear filtering

texts:
I am wondering if a methodology(ies) exist for essentially 'filtering' out a tone or interference, WITHOUT resorting to classical linear filtering techniques. (band pass filter, etc for example).

The reason for this question is because inherent in any classical linear filtering operation, I am at the end of the day linearly combining samples among themselves, and while this might do well to remove a tone of another frequency that I am trying to remove, it also has the effect of 'smearing' my signal in time. This smearing effect causes information such as the 'sharpness' of the signal to be lost, and in this specific case, the 'sharpness' carries information such as start of signal energy, etc, that I use later.

What method(s) exist (do they exist?) for this goal?

I should also add, please assume omnidirectional case, (I only have one sensor), but if assuming a case with multiple sensors lends itself to a solution, feel free to use that assumption as well.

Thanks

EDIT#1:

I must say that for the tonal case,  I think it might be straight forward as many of you have indicated, vis-a-vis estimating the tonal parameters, reconstructing, and subtracting. The actual problem though, is a little more complicated...

The unwanted signals can range from anything like a group of whales having a 'conversation' to a boat's motor overhead. In such cases no known robust model exists, all I know is that my signal of interest is among all those other energies - they are seperable in frequency, but I need a way to remove them as well, without filtering because that removes my sharpness.

FYI, I think this falls into the realm of non-linear filtering or other techniques from different areas? I should also add that I have a lot of computing power.


texts: Is sub-sampling delay estimation possible without interpolation?

texts: I want to estimate time delay between two complex narrow band signals. I can use any correlation method. There are several methods in the literature about sub-sample time delay estimation. Generally this subs-sample delay estimation is done through phase estimation. My question is: are these sub-sampling methods using phase requires interpolation or not? I was expecting that we are interpolating in Fourier domain and using phase information to estimate time delay. Any comment would be appreciated.


texts: Is it wise to use mathematical morphology and Gabor filter to extract blood vessels in the same algorithm?

texts:
I'm just curious, is it really wise to use mathematical morphology and Gabor filtering in the same code to extract blood vessels?

Like in this case, I use morphology to get the vessels and using Gabor filter to further enhancing brightness of my vessel.


texts: How to solve this question on communication lines?

texts:
The Question asked is as follows:
A communication link operates at a bitrate of 20 Mbps. The receiver is expected to perform at a $BER$ of $10^{-6}$ or better using $QPSK$ modulation scheme.
[Note: $BER$ for $QPSK$ is given by $Q(\sqrt{2Eb/
N_0})$. Use the approximation $Q(x)≈ \frac{1}{2} e^{\frac{{−x}^2}{2}}$]
(a) Determine the minimum $\frac{Eb}{N0}$ in dB needed to achieve the desired $BER$ performance.
(b) Given that the noise figure of the receiver is 8 dB, determine the receiver sensitivity
in dBm.
(c) If the required link margin is 20 dB and the transmit power $P_{Tx} = 0.5 W$, determine
the maximum possible path-loss, $L$, in dB
(d) Using a carrier frequency of 5 GHz and assuming free-space propagation model,
determine the maximum possible range in $m$ for this system.
I attempted to solve the question, but my answers didn't match with the given answers. Any help on how to go about solving this would be highly appreciated. Thanks.


texts: Difference between texture and illumination in image

texts:
I am new in image processing.

What are texture and illumination in image?
How can we define these terms?


texts:
I am looking for algorithm(s) to solve the following problem:
Given a noisy .wav sound capture (some wind + friction noise on the microphone), how to detect the BPM of a soft drum beat?

I have attempted googling the subject, but the results are quite poor, due to high amount of mp3 related software for both analysis and fingerprint id generation. None of them supply information about how to actually do it.

I am aware of algorithms to remove the noise, but that still leaves me with the problem of detecting BPM. And depending on how the BPM problem is solved, it's possible that I don't even need to denoise (since drum tends to be in the lower frequencies and noise in higher, a simple low-pass might be sufficient pre-processing).


texts: Matlab's invfreqs won't fit at low frequencies

texts:
I'm wondering of anyone can explain why invfreqs() is unable to fit a polynomial to the data in the image below.  The red line is the measured frequency response of an analog system.  I should mention that this data is a logarithmic sweep from 100 Hz to 100 GHz. The blue line is the 11th order invfreqs fit (11th order in both numerator and denominator).  11th order does the best job, but as you can see, there is a massive offset from DC to a few MHz.  Everything from 1 GHz and above is a stellar fit.

Any idea why the low freqs don't fit?

UPDATE:  Problem solved. The frequency vector must be normalized to keep the computation accurate.  I divided it by 1e6.  Then the b,a coefficients are converted into zpk form.  Multiply z and p by 1e6.

%% estimate transfer function using invfreqs()

H = freq_resp;
samples = length(H);

w_norm = w/1e6;

wts = 0:samples-1;
wts = exp(-wts);

[b,a] = invfreqs(H,w_norm,15,15,wts,50);

[z,p,k] = tf2zp(b,a);

zs = z*1e6 ;
ps = p*1e6 ;

[bs,as] = zp2tf(zs,ps,k);

sys_pred = tf(bs,as);
figure
bode(sys_true,'r',sys_pred,'b')


texts: Ramp function as derivative in frequency domain?

texts:
It is said that to get Laplacian of Gaussian in frequency domain, we may multiply the Fourier transform of Gaussian with two differentiating ramp function (1 ramp gives 1 order of derivative).

The description from the material that I was following:

And the file can be found here

So how could we possibly get derivatives by multiplying ramp functions in the frequency domain? How does the math work here?


texts:
What do you call it when you add zeros between the elements of a vector?

let's say you have x =[1 2 3 4]

and x_2 = [0 1 0 2 0 3 0 4]

what do you call this process?


texts: Measuring the similarity of two waveforms

texts:
I have two approximately sine truncated waveforms of different amplitudes and I want a numerical measure of how similar they are to each other. I am sending a sonic pulse train whose multiple reflections I am receiving. I want to know which of those reflections most closely correlate with each other in terms of wave shape.
Any standard methods?


texts: Audio Watermarking - What is typically embedded?

texts:
Just starting research in audio watermarking and have read many papers regarding audio watermarking but I am still a little confused on what exactly gets embedded into the signal. For example, one paper I read embedded an image or written signature into the audio.

My question is what typically gets embedded in the signal using audio watermarking? Is it just random data, a picture, someone's name? Thank you!


texts: FMCW SAR chirp return from target modelling

texts:
I have to Model a chirp signal that was return from a Distance D from target and then multiply with the reference transmitted signal. Which will give me beat frequency signal and then fft to get the range.
For this I was searching around and find a pulse LFM SAR implementation.

Nf=2^10; %number of points for freq. range
fc=242.4*1E6; %stopfreq
B=133.5*1E6; %bandwidth
freq=linspace(fc-B,fc,Nf);

function [resp] = my_sinthSISO_planresp( pos,freq,tar,tar_amp )
%MY_SINTHSISO_PLANRESP calculates synthetic SISO at position pos response
%in frequency range freq on target tar(target postion), tar_amp: RCS like.
%Targets and Array must be in the same (image) plane, which XY-coordinate
%are represented as complex values.

if nargin < 4, tar_amp=ones(size(tar)); end % tar_amp = 1 if not given

%speed of light
co=299792458;
%redistribution of input data
tar_dim=numel(tar);
tar=reshape(tar,[1,1,tar_dim]);
tar_amp=reshape(tar_amp,[1,1,tar_dim]);

freq_dim=numel(freq);
freq=reshape(freq,[1,freq_dim,1]);

pos_dim=numel(pos);
pos=reshape(pos,[pos_dim,1,1]);

%% algo
%targets-elements radial-distances matrix
D=abs(repmat(tar,[pos_dim,1,1])-repmat(pos,[1,1,tar_dim]));
D=repmat(D,[1,freq_dim,1]);
%wavenumber
k=2*pi*repmat(freq,[pos_dim,1,tar_dim])/co;
%target amplitude
tar_amp=repmat(tar_amp,[pos_dim,freq_dim,1]);
%SISO response
resp=tar_amp.*exp(-1j*k*2.*D);
resp=sum(resp,3);

end

Question 1: Is resp in a frequency domain? since its dependent on k and k is frequency related.
Question 2: if its freq domain then why the inverse fft does not give the desired chirp signal. while its showing:
Question 3: If it doesnot represent chirp in frequency domain. Then what it represents?

Question 4: How can I use this data set to model FMCW. Which means multiplication of the chirp at reference and after shift?

I know its quite alot and very basics but I got confused. Thanks.


texts: How to verify if my cross-correlation algorithm works?

texts:
So I have a signal generator which can produce two separate signals.

I use a spectrum analyzer to read these signals, and can interface with Python.

I wrote a Python code which takes in these two signals, perform the FFT of both, take the complex conjugate of one of them, and multiply together which gives me the cross-correlated frequency spectrum?

My problem is, I want some validity checks to see if this works or not. The main objective is noise floor reduction, but I am not sure how much the noise should drop by, etc. This signal generator is connected to the spectrum analyzer using BNC coaxial cables and that is it. I am sending in simple sinusoidal waves of frequency 2kHz.


texts: I still get aliasing of my signal, even though I stay below the Nyquist frequency. What am I missing?

texts:
Thanks in advance for any help with my newbie question. I've been consulting quite a large number of resources in order to understand how I can create and playback a nice sawtooth without aliasing, but I feel something is really not clicking, since explanations I've been reading have only been confusing me further.

It's my understanding that aliasing should not occur if my signal does not contain any sinusoidal components above the Nyquist frequency. So, to be safe, I create a sawtooth using the following code:

            fs = 48000;
            a = linspace(0,(fs/20),fs);
            SinAudio = sin(2*pi*a+0.5*pi);

            for i = 1:(length(amaxi)-1)
                SawElement{i} = linspace(-1, 1, (amaxi(i+1)-amaxi(i)));
            end

            % Make first sawtooth element
            Saw = linspace(0,1,amaxi(1));

            % Concatenate rest of elements
            for i = 1:length(SawElement)
                Saw = [Saw, 0, SawElement{i}];
            end

In this case, what I perceived as aliasing is apparent as a lot of clearly audible low frequencies, that sort of sound bubbly, especially when I sweep the signal.

So, the base frequency is way below my sampling frequency, and the way I concatenate the elements of the sawtooth, I take two samples to go from 1 to -1. Intuitively, I feel like this is where it goes wrong, because if you were to create a signal like that from individual sinusoids, you would need A LOT of harmonics, which would be why I get aliasing..

HOWEVER, there should be no more harmonics there than when I create a sine wave like this:

fs = 48000;
a = linspace(0,(fs/4),fs);
x = sin(2*pi*a+0.5*pi);

Which is a sinewave with a frequency of half the nyquist frequency, and therefor I wouldn't expect aliasing... I think?

So, my questions are:
1) Why does a generated signal with maximum frequency at half the nyquist rate still alias?
2) How do I create a bandlimited sawtooth or squarewave? I tried using wavetable synthesis, and although I haven't tested the resulting signal yet, I don't see how bandlimited wavetable synthesis would give me any different results from the two pieces of code I described above; ultimately, there's going to be a point where my signal goes from 1 to -1 over the course of 2 or 3 samples, no matter how I generate my signal..

Thanks in advance, I've really been breaking my head over this....

Paul


texts: Projecting a vector to another to detect the sparse values of such vector

texts:
Assuming we have sparse vector of length $N$ such as $X = [0,1,0,-1,0,1,1,0]$ which has some non-zeros values. The vector $x = iFFT(X)$ is convoluted with another vector $h$ resluting $y = h*x$. Suppose that $h$ and $y$ are known, is it possible to recover the sparse vector $X$ using such methods of sparse vector estimation, such as compressive sensing or any other method?

$NP$:
I have suggested to take $FFT$ for $y$ resulting $Y = FFT(y) = HoX$, where $o$ is the element-wise multiplication between $H$ and $X$. Then based on $Y$ I think it's easier to recover sparse vector $X$. But the issue when using element-wise division, I can't get the position of each non-zeros element in $X$. On the other hand, I don't know how to implement compressive sensing to recover $X$.

Thank you in advance


texts: Digital demodulator (SDR) for receiving legacy analog modulated wave

texts:
I understand that in SDR, after the ADC the signal processing is fundamentally based on I/Q modulated waves where any discrete value that can be transmitted corresponds to a symbol which is represented as amplitude and phase through I and Q signals.
How can a modern SDR receive legacy AM waves?
The SDR uses I,Q components whereas legacy AM doesn't have any such component applicable.
I presume SDRs are able to receive and demodulate existing legacy AM,FM transmissions - please correct/add.
Do current SDRs receiving AM waves switch to a different technique wherein it doesn't use I/Q demodulation?


texts:
I'd like to validate changes to our app music player in our Continuous Integration unit tests (more like integration tests, really). My thoughts were to start with an audio file that contains only a pure tone, or a series of pure tones. Our automated test would stream this file over the internet and play it through our player. We would capture the output using MTAudioProcessingTap. This output would go through and FFT and the pass/fail test would be check for the existence of only the original frequencies.

Would this work? I'm thinking this would catch pops and glitches, but not necessarily discontinuities in the output.


texts:
Suppose you have the ability to inject any arbitrary waveform into a piece of analog rf hardware and collect and digitize the output for analysis.

If you wanted to characterize/estimate the transfer function of the hardware, what waveforms would you use?


texts: FIR Filter design, using different filter length than signal length

texts:
Given some neurophysiological application I am filtering in real time data of limited length (e.g. in epochs of 100 ms) using a sampling rate of e.g. 100 Hz.
Due to this short time segment I am trying to implement a FIR filter. This filter should be of high order, and therefore have a length of more than 100 samples.
This would lead to the fact that the "to be filtered signal" has a different length than the estimated filter itself. Using the scipy.convolve function I can then specify different methods how to convolve the signal (full, valid, same).
My question is quite trivial, and related to lack of experience in the field, but does it make sense at all to have different filter and signal length?
Is there literature how to interpret then the filtered signal? Since the convolution needs to define some padding for the signal itself the filtered signal should be only in a certain range interpretable right?
EDIT:
Hereby I want to present an example. I plot a random signal in orange, and a FIR filter with a filter length of 10 s (estimated using mne.create_filter):

My understanding problem is visualized in this context: the signal itself is much shorter than the filter. This leads obviously to the fact that the filtered signal has a longer length. How can I now interprete the filtered signal with respect to the timing of the original signal that had fs=128 Hz? Which time points of the filtered signal correspond to the original signal?


texts: Find the $\mathcal Z$-transform of this function?

texts:
I need to find the $\mathcal Z$-transform of $x(n)=a^{-n} u(n)$. Assume, $a$ is a positive constant , but the power of $a$ is negative.

Looking at the transform table, I found that $\mathcal Z$-transform for $x(n)=a^{-n} u(n)$ is not there, so I tried to found it by on my own, but I am not sure if my answer is right, please can anyone check it and correct me, in case?


texts:
Here I understand that the first term is a simple convolution of the input signal with first volterra kernel that I have acquired through farina sweep method, my question is if the second term is the convolution of second volterra kernel with our input signal raised to the power of two? Would the code using matlab, of the first and second terms be correctly written as `

outputsignal1=conv(voltkernel1,inputsignal.^1); %basic convolution with system impulse
outputsignal2=conv(voltkernel2,inputsignal.^2); %convolution of second kernel with squared signal
summedsignal=outputsignal1+outputsignal2;  %dynamic nonlinear output signal with memory


texts: phase shift between two images (matrices)

texts:
I'm trying to find the phase shift between two 100x100 matrices, s and p. To start with I'm using two sine waves (2D) slightly shifted. I manage to compute the correct amplitude spectrum (using the Fourier Transform of the cross correlation between s and p) but the phase value at the peaks is not the one that I have set. I'm not using cpsd because I did not manage to understand the correct way to use it for 2D matrices. The manual is very confusing and I couldn't get any help online.

When I do the same procedure for 1D sine waves, using cpsd or fft of xcorr employing the welch method, I get both the amplitude and phase right. I suppose that conceptually I'm doing the right thing, but there's sth about how to get the phase shift between two matrices that is bugging me.

I really would appreciate any tips on how to compute the phase shift between matrices. Some real piece of code would be lovely.

thank you in advance for your time


texts: Oversample or undersample correlated noise?

texts:
I wish to measure nontrivial (non-Gaussian) broadband current noise for its distribution and power spectral density.
The noise is amplified (with transimpedance amplifier) with some characteristic risetime (tau ~30 usec), that I assume makes the signal correlated with the same timescale. The noise itself may be correlated with different timescale (~1 usec to 1 msec).

Is it best to oversample (sample frequency > 1/tau) or undersample (< 1/tau) the signal?
Does it really matter and what is the difference/effect?

I thought that slight undersampling draws directly from the distribution and the distribution is not biased to "double counting" (measure distortion) of the signal, however a colleague advised to overasample, thus no information gets lost.

To make it clear, the signal is random noise, and the purpose is to properly characterize it (so I'm not sure how Nyquist theorem fits in or relevant to it. Maybe other sampling theorems may be useful).


texts:
For evaluating the performance of algorithms that use double-precision floating point numbers, or to pre-calculate double-precision floating point data to the best precision, it would sometimes be useful to evaluate the discrete Fourier transform (DFT) and the inverse discrete Fourier transform (IDFT) at arbitrary precision greater than the precision of double-precision floating point numbers. Python's mpmath library would be a well-suited framework on which to build an arbitrary precision fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT) library.
Is there such an implementation available? If no, then this is a request to develop and present a reference implementation.
(Hint:  chirp-z transform might be the easiest way to deal with any DFT length that is not a power of two.)


texts: Perfect reconstruction ( Single Channel )

texts:
I was trying to solve the following problem .

  Given an input x[n], consider upsampling by 2, followed
  by interpolation with a filter having z-transform H(z) for magnification of the signal.
  Then, to recover the original signal size, apply filtering by a decimation filter G(z) followed
  by downsampling by 2, in order to obtain a reconstruction ˆx[n].

using the results for upsampling as $Y(z) = X(z^2)$

and downsampling as $Y(z) = X(z^{1/2})+X(-z^{1/2})$

I obtained

$\hat{X}(z) = 1/2\bigg(G(z^{1/2})H(z^{1/2})X(z) + G(-z^{1/2})H(-z^{1/2})X(-z)\bigg)$

however, the solution is listed as ,

$\hat{X}(z) = 1/2X(z)\bigg(G(z^{1/2})H(z^{1/2}) + G(-z^{1/2})H(-z^{1/2})\bigg)$

Where have i gone wrong ?

I solved in the following process ..

$A(z) = X(z^2)$

$B(z) = H(z)A(z)$

$C(z) = G(z)B(z))$

$\hat{X}(z) = 1/2[C(z^{1/2})+C(-z^{1/2})]$

and substituted them in a chain.


texts: Why does this Moiré pattern look like this?

texts:
I was making some gifs of Mobius transformations in Matlab, and some strange patterns began to appear. I'm not sure if a deeper knowledge of the filetype/algorithm is needed to understand this phenomenon, but I thought that there could perhaps be a purely mathematical explanation. The image is obtained by coloring the complex plane like a checkerboard, and then inverting it by taking the reciprocal of the complex conjugate. Here is the math psuedocode for the image with a given zoom $k$:

\begin{align}
&\mbox{checkerboard}:\mathbb C \to\{\mbox{black},\mbox{white}\}\\
&\mbox{checkerboard}(z):=\begin{cases}
\mbox{black} & \mbox{if }\lfloor\Im(z)\rfloor+\lfloor\Re(z)\rfloor\equiv 0\mod 2\\
\mbox{white} & \mbox{if }\lfloor\Im(z)\rfloor+\lfloor\Re(z)\rfloor\equiv 1\mod 2
\end{cases}\\
&\mbox{image} = \{z\in\mathbb C:|\Re(z)|,|\Im(z)|\leq 1\}\\
&\mbox{color}:\mbox{image}\to\{\mbox{black},\mbox{white}\}\\
&\mbox{color}(z):=\mbox{checkerboard}(k/\overline{z})
\end{align}

And here are the pictures for $k=1$, $k=50$, and $k=200$. The resolution of each picture is $1000\times 1000$. I have no background in signal processing, but I would be eager to learn things!

EDIT:

More specifically, why does the Moiré pattern 'sync up' with the resolution of the picture at certain points?
Can the Moiré pattern be predicted?


texts: Definition of Ideal Low pass filter (Time Continous)

texts:
Actually I got confused about definition of frequency response of Ideal low-pass filters because in some books they mention that phase of H(f) should be linear in pass band and it's value is $ -2\pi t_{0}$

However in other books they mention that phase should be zero in pass band.

So I want to know their differences and where does $ -2\pi t_{0}$ comes from?


texts:
My question is related to this article: https://www.dsprelated.com/showarticle/192.php

I think I understand mostly everything until this sentence:

"The directions in which the impulses are pointing show the relative phases of the spectral components".

Please, note that I understand (and tested in Matlab) that the DFT of a pure sine wave of frequency f has a negative peak at the positive freq f and a positive peak at freq -f in the imaginary part Im(fft(x)) and no peaks in Re(fft(x)). Cosine has only positive peaks in Re(.) and nothing in Im(.) etc.

However, I have trouble understanding this story about phases, can someone help me understand this sentence:

"The directions in which the impulses are pointing show the relative phases of the spectral components".


texts: Detecting inset images/videos from a video

texts:
I want to detect multiple rectangles in a video, each of which features a different content, but running simultaneously. This is typically seen during debates on news channels. For example, in the image below, there are 5 rectangles (each containing a person) that are of interest to me and I want to keep a track of them. I am planning to use C++ and OpenCV. Any help will be appreciated.
Thanks.


texts: Single-sided Z transform with difference equations and the system function

texts:
I am working on this problem:
Given an impulse response,

find the system function,
find the difference equation representation,
find pole-zero plot,
find output $y[n]$ if the input is $x[n] = 0.25^n * u[n]$

Here is what I have so far (MATLAB code)

Find H(z)
$$H(z) = \frac{5}{1 - 0.25z^{-1}}\quad,\qquad|z| > 0.25$$

For part 2 I used the fact that
$$H(z)=\frac{Y(z)}{X(z)}$$
I cross multiplied, then took the inverse $\mathcal Z$ transform and got:
$$y[n] - 0.25y[n] = 5x[n]$$

For part 3 I used the zplane function:

 b = [5 0];
 a = [1 -0.25];
 figure; zplane(b,a);

Part 4 is where I am getting confused. I took the one-sided $\mathcal Z$-transform of the difference equation to get
$$Y^+(z)-0.25[y[-1]+z^{-1}Y^+(z)] = 5X^+(z)$$
Assuming the system is casual (due to unit-step function in impulse response and input $x[n]$, also no initial conditions given in the problem),
$$y[-1] = 0$$
So, after rearranging, $$Y^+(z)=\frac{5}{1-0.5z^{-1}+0.0625z^{-2}}\quad, \qquad |z| > 0.25$$
Using the residuez function for partial fraction decomp

 b1 = [5 0 0];
 a1 = [1 -0.5 0.0625];
 [R p C] = residuez(b1,a1)

 R =

      0
      5

 p =

     0.2500
     0.2500

 C =

      0

Then taking the inverse $\mathcal Z$-transform of the terms generated by residuez,
$$y[n] = 5(1/4)^nu[n]$$
However, when I check, this is not the case. For the check, I first create a MATLAB generated $y[n]$ sequence using filter

%MATLAB check for part 4
n = 0:50;  % compare first 50 samples
x = (1/4).^n;
y = filter(b,a,x);
Then, create my check sequence calculated from the inverse $\mathcal Z$-transform and the residuez function

ycheck =  5*(0.25).^n;
error = max(abs(y-ycheck))
The sequences are not the same, and error nowhere close to nominal. Can anyone tell me where I went wrong please?


texts: Is there any tool to anlyze audio specturm of a song and give the average visually ?

texts:
I have a lot of songs I made over the 20 years and was thinking it would be interesting to see how the overall EQ changed on the various songs I made.

Any tool out there to help me with something like this ?