texts: Image processing: using PCA for 3d object recognition

texts:
I read this lecture on using PCA/eigenfaces: http://www.cse.psu.edu/~rcollins/CSE486/lecture32.pdf

And I want to use this general procedure to recognize a 3D clay model that I took pictures of from multiple angles.

I'm having trouble conceptualizing how it's actually implemented in software. I've "trained" the initial mean and eigenvectors for using ~40 images offline. So I have new input images that will have my 3D model in it (or wont).

My confusion: I'm worried that there's so much else going on in my input image that applying the projection + thresholding equation isn't going to work. It almost seems like I need to sub-image my input many times and try over and over again to get the right projection (almost like multi-template matching). Or maybe I can detect blobs first, create sub-images from the blobs, then try the projection.

Because I have multiple angles of the model I suppose that handles rotational invariance, but what about scale invariance? Seems like whatever subimage I get from my input, I have to try at different scales. Or I'm not understanding something fundamental about how to use PCA.

Thank you for any help; for reference I would like to implement this in OpenCV (Java-based).


texts: Confidence Interval For Frequency of Peak in Averaged Periodogram?

texts:
I am averaging periodograms from sets of ~15-25 experimentally measured timeseries all of the same length.  To do this I am just summing the power in each frequency bin over all periodograms and then dividing by the total power over all the periodograms (I believe this is similar to Welch's Method but here I have many timeseries instead of a single long timeseries).

In many of these averaged periodograms I see a distinct peak that is called as statistically significant by Fisher's G-test.  What is a good way to get a confidence interval or other interval estimate (e.g. standard-error of the mean) for the frequency of this peak?  Is there some method that takes into account both sample size and width of the peak (maybe width at half-max)?

Add note below to provide more detail:

These are real timeseries data from measurements of gene expression in cells (each timeseries is the gene expression in an individual cell sampled every 6 minutes for ~10 hrs).  I have been averaging the individual periodograms  to create an averaged periodogram (like Bartlett's method) and want to be able to determine two things:

i.) Is a power at a frequency in the averaged periodogram statistically significant (i.e. unlikely to be due to chance)? This is what I am currently using Fisher's G-Test for but if there are better methods I am open to trying those (provided I can understand them and they are relatively easy to use).

ii.)  What is an interval estimate for the frequency of a statistically significant power in the averaged periodogram?  The bin with the greatest power in these averaged periodograms is typically a "peak" (when a line is plotted connecting the power values across bins) with neighboring bins also showing some increased power relative to other parts of the periodogram.


texts: length of window and overlap rate in STFT

texts:
I want to use STFT to analyze my signal and am wondering what are differences between two solutions:

Use short windows (for ex. 256 samples window)
Use longer windows (to get higher resolution in freq. domain) but let the windows overlap more (For ex. I use 1024 samples window, but let 1023 samples overlap)

Thank you :)


texts: In Convolutional Neural Nets, what do convolutions look like?

texts:
The early stages of the Convolutional Neural Networks are performing classical convolutions with a certain kernel size on the input image. Is is possible to express in common terms the type of convolutions they are generally performing (lowpass, highpass, differential, anisotropic...), or are they completely application-dependent/impossible to interpret ?

Alternatively, can you show some sample kernel ? I would be very grateful.


texts: Find a small image pattern in another similar image

texts: I have this pattern in the first image which I want to find in the second image. However the two images are totally different so I just need to find the closest set of indices where these two images match. So how to go about it? Are there existing algorithms doing the same?


texts: Amplitude of frequency bin in FFT doesn't match time-domain amplitude

texts:
I have a sinusoidal current that I am sampling at about 357k SPS. The current signal is about 3A pk-pk @ 750hz (top figure). I am wondering why when I take the fft using numpy/scipy's fft function, the amplitude of the fundamental frequency doesn't match.

I have a sample size of 8000 so I divided the fft by 8000 and multiplied by 2 to get the results below.

Current and fft calculation:

# current sense
df['current'] = (df['Viout'] - 2.5)/20/0.003
df['current_ft'] = np.fft.fft(df['current']) / 8000 * 2

Plot:

plt.subplot(211)
plt.title('current')
plt.xlabel('time [s]')
plt.ylabel('current [A]')
plt.plot(df.t_1, df.current)

plt.subplot(212)
plt.title('current fft')
plt.xlabel('frequency [Hz]')
plt.ylabel('|F(f)| / 8000 * 2')
plt.plot(df.freqz, np.abs(df.current_ft), 'o')
plt.xlim(-ax, ax)

plt.tight_layout()
plt.show()


texts: Gnuradio fifo cntrl timed out looking for acks

texts:
I am using USRP N210 SBX for my experiment. It is a simple TX RX system with a tuning center frequency. I have uploaded the flowgraph.
I am connecting to the USRP using a remote server
CPU MHZ : 2300.045
cache size : 512KB
I get the error " fifo ctrl timed out looking for acks" sometimes. Sometimes it executes. I have copied the whole error
linux; GNU C++ version 5.3.1 20151219; Boost_105800; UHD_003.009.002-0-unknown

Xlib:  extension "XInputExtension" missing on display ":1".
-- Opening a USRP2/N-Series device...
-- Current recv frame size: 1472 bytes
-- Current send frame size: 1472 bytes
Traceback (most recent call last):
  File "/home/banerjee.152/rf_spectrum.py", line 179, in <module>
    main()
  File "/home/banerjee.152/rf_spectrum.py", line 167, in main
    tb = top_block_cls()
  File "/home/banerjee.152/rf_spectrum.py", line 76, in __init__
    channels=range(1),
  File "/usr/lib/python2.7/dist-packages/gnuradio/uhd/__init__.py", line 122, in constructor_interceptor
    return old_constructor(*args)
  File "/usr/lib/python2.7/dist-packages/gnuradio/uhd/uhd_swig.py", line 2671, in make
    return _uhd_swig.usrp_source_make(*args)
RuntimeError: RuntimeError: fifo ctrl timed out looking for acks

I have read the previous threads on this topic and increased the buffer to max
sudo sysctl -w net.core.rmem_max=50000000
sudo sysctl -w net.core.wmem_max=50000000
However I still find the same issue sometimes.
I repeated the whole experiment in my local host laptop. I dont find any error. The major difference is
cpu MHz     : 2840.578
cache size  : 8192 KB
Is the cache size the real issue or can I fix this problem while working with the remote server?
Any understanding will be greatly appreciated


texts: Using the system function $H(z)$ to the solve the output of a DT LTI system

texts:
I have a question which asks me to find the magnitude and phase response of

$$
H(z) = \frac{1}{1-0.5z^{-1}}
$$
for which I solved by doing
$$
H(e^{j\Omega}) = \frac{1}{1-0.5e^{-j\Omega}} = \frac{1}{1-0.5\cos(\Omega)+j0.5\sin(\Omega)}
$$
$$
|{H(e^{j\Omega})}| = \frac{1}{\sqrt{(1-0.5\cos(\Omega))^{2}+(0.5\sin(\Omega))^{2}}} = \frac{1}{\sqrt{\frac{5}{4} - \cos(\Omega)}}
$$
$$
\angle H(e^{j\Omega}) = -\tan^{-1} \left(\frac{\frac{1}{2}\sin(\Omega)}{1-\frac{1}{2}\cos(\Omega)}\right)
$$
But, I am unsure how to solve the second part of the problem, which asks me to compute the system response to
$$
f[k] = A\cos\left(\frac{\pi}{2}k +\frac{\pi}{10}\right)
$$


texts:
Suppose that a system induces a phase of -90 degree, how do I practically measure the delay in the signal coming out of the system? I think time lag would be a more useful measure than phase.

Clearly these have different units, what is the conversion factor that tells me how much time lag is induced due to the system?


texts: Decreasing Sample Rate in DFT (FFT) for Audio Analysis

texts:
I'm working with audio data analysis through FFT algorithm.
My example audio is a sine wave at 440 Hz and 44100 Hz sampling rate.

FFT methods in programming, like scipy.fft(y), require a vector of samples and a number of points/samples to analyze.
Sample rate is used to determine the width of each frequency range/bin in the resulting data vector (FS/N).

I've read that to increase frequency resolution of FFT results one should dicrease sampling rate and increase window size (number of samples).

I understand that I can change the N value (samples count) that I pass to the algorithm.

But if the audio is at 44100 Hz how can I dicrease it?
Would the resulting frequency bands reflect incorrect information about magnitudes of the FFT data?

Or should I change the input audio data somehow (manipulating directly the samples in the vector) before FFT so the data is passed with lower rate (say at 22050 Hz)?

Thanks.


texts:
I use the following formula for DoG:

$$\frac{1}{\sigma} (\frac{x^2}{2\sigma^2}-1.0) e^{\frac{-x^2}{2\sigma^2}}$$

What is the relationship between this formula and the difference of two Gaussian functions?

Can someone show me how to apply DoG on a 1D signal in MATLAB? My input can be considered as a 1D real-valued function, and I applied the above formula to detect blobs. The results make sense except that I don't know why I get many local minima in a relatively flat region (i.e., function values at points in the region almost stay constant). Since I don't have enough signal processing/image processing background, I want to see what a constant signal look like after being convolved with a DoG filter.


texts:
From this Wikipedia link, the fastest fequency for digital frequency is 0.5cycle/sample. It can be rewritten as 1cycle/2sample.

What does it mean?
Does it mean that at least two pixles are required to form one cycle?
Is there any requirement on the values of two samples?
One must be 0 and the other must be 255? If so, why?

Edit:
The reason I have such question is that different cutoff frequencies result in different blurring degree for filtering in frequency domain.

From the link:http://www.imagemagick.org/Usage/fourier/#fft_partial

The figure shows the comparison. I hope to estimate blurring degree from cutoff frequency in digital frequency. Maybe I can estimate corresponding standard deviation of Gaussian filter from cutoff frequency.


texts: Nonlinear time-invariant frequency doubler

texts:
Does there exist a causal nonlinear, time-invariant system that maps each input function $k \cdot\cos(fx)$ to $k \cdot\cos(2fx)$ for all choices of $k, f \in \Bbb R$?
If so, can one represent this system via some sort of differential equation, or in discrete terms, as a nonlinear difference equation?

This would be similar to the second Chebyshev polynomial, but that only works for $k=1$.


texts: Regarding custom wavelets and their validity

texts:
I've recently started working with wavelets, and want to focus primarily on Discrete Wavelet Transform. For experimental purposes, I use the PyWavelets library in Python. In one of my tests, I tried the following code.
Basically, I randomly sample the coefficients for the lowpass and highpass filters from a normal distribution, and use those coefficients to plot the scaling and wavelet functions.
import matplotlib.pyplot as plt
import pywt
import numpy as np

fig = plt.figure()
fig.subplots_adjust(hspace=0.2, wspace=0.2, bottom=.02, left=.06,
                        right=.97, top=.94)

dec_lo, dec_hi, rec_lo, rec_hi = np.random.normal(0.0, 0.5, (4, 2)).tolist()
print(dec_lo, ' ', dec_hi, ' ', rec_lo, ' ', rec_hi)
filter_bank = [dec_lo, dec_hi, rec_lo, rec_hi]
wavelet = pywt.Wavelet(name="myWavelet", filter_bank=filter_bank)

phi, psi, _, _, x = wavelet.wavefun(level=5)

color = 'red'
ax = fig.add_subplot(1, 2, 1)
ax.set_title(wavelet.name + " phi")
ax.plot(x, phi, color)
ax.set_xlim(min(x), max(x))

ax = fig.add_subplot(1, 2, 2)
ax.set_title(wavelet.name + " psi")
ax.plot(x, psi, color)
ax.set_xlim(min(x), max(x))

plt.show()

A run of the above code gave me the following plot (Since values are randomly sampled, different runs will give different plots):

My question is, can it be shown that these plots are valid scaling and wavelet functions? Or are they simply meaningless plots generated for the sampled coefficients?


texts: How does spectrum analyzer analyze high frequency?

texts:
I have a question that is regarding spectrum analyzer. According to nyquist theorem, sampling rate  should be twice of maximum frequency of signal to reconstruct again. If there is a spectrum analyzer that is range 2kHz to 20GHz, is sampling rate equal 40GHz?
How to analyze high frequency component? does it pull IF frequency?


texts: LTI system output given input and frequency response

texts:
The question I'm trying to understand is as follows: A linear time-invariant continuous-time system has the frequency response function $$H(\omega)=\frac{1}{j\omega+1} $$

Compute the output response $y(t)$ for $-\infty<t<\infty$ when the input $$x(t)=\cos(t),\quad -\infty<t<\infty$$

I know the answer, but I don't understand how to get there.

I know $H(\omega)=\frac{Y(\omega)}{X(\omega)}$, and from that $y(t)=h(t)*x(t)$. Taking the inverse Fourier transform of $H(\omega)$ I get $e^{-t}$. If I try to compute the convolution I just get a divergent integral. Is this wrong? The answer is supposedly $$y(t)=|H(1)|\cos(t+\angle H(1))$$ but this seems to come from nowhere. Can someone please explain?


texts:
Can you advise me on good inhomogenity measurements in an image?
I am trying to find the borders where homogenous zones touch inhomogenous ones in a image.

I am going to work with sick lung ct images like this:

I am looking for different ways to measure inhomogenity in the texture in the image.

Thank you


texts:
Sorry if I ask a basic question but I am a bit lost.

I have a data set in frequential domain, I need to go to time domain to perform time-gating and then going back to frequencies. In order to do this, I compute an inverse Discrete Fourier Transform.

Shannon theorem tells us that two times the difference between the greatest and the least frequency, respectively noted as $f_N$ and $f_0 ( \ne 0)$, contained in the signal must be less then the sample-rate in time domain ($F_t$). Do I need to do some zero padding to avoid aliasing ?

Without zero padding on frequential domain, what is the relation between $F_t$ and $f_N, f_0$ ? My guess is that $F_t=f_N-f_0$. Is that correct ?

If I must do some zero padding, should I do it this way :

$$(0,\dots,0,f_0,\dots,f_N,0 \dots,0)$$

With enough zeros such that the first frequency contained in the padded signal is $0$ and the last is $2 f_N$ ?
If my guess above is correct, then the new sample rate is high enough to avoid aliasing, because it will be $F_t=2f_N>2(f_N-f_0)$.

Thanks for your help.


texts: What is the relationship between a digital audio signal and sound pressure level?

texts:
I'm playing a pure tone on an iPad through a pair of "standard" earphones.  From a programming perspective, the output level is in the range of $0$ to $1$ where $1$ means "full volume".

I've observed using a sound pressure level meter that when playing a pure tone at $1000 \text{Hz}$ with output level in $0$ to $1$ with step size $0.05$, that the relationship between output level and sound pressure level is logarithmic.

However, I'm not sure if this is correct or an issue with my use of a sound field sound pressure level meter (instead of a coupler).

Should doubling the output level double the sound pressure level?  That is, should the relationship be linear?  It makes intuitive sense to me that as the output level or voltage is increased, the earphone speaker will vibrate more, and thus the sound pressure will also increase.


texts:
Three methods of estimating maximum signal level at a node:

Worst Case Min-Max (L1-Norm)
Mean Square (L2-Norm)
Frequency Domain (Linfinity-Norm)

In matlab i found the code for L1, L2

Hfull=[h(r,:),conj(fliplr(h(r,2:end-1)))];
L1(r)=sum(abs(ifft(Hfull)));
L2(r)=sqrt(mean(Hfull.*conj(Hfull)));

For Linf is correct the below?:

Linf(r)=max(abs(Hfull));


texts: how were DTMF frequencies determined exactly and how can I extend them?

texts:
I've been playing with DTMF frequencies and the Goertzel algorithm, and would like to know exactly HOW the DTMF frequencies were determined. I understand that they were chosen such that none of the frequencies would be equal to sum and difference frequencies, or intermodulation products, but I can't find any information on whether the frequencies were simply determined through experimentation, mathematical calculation, etc… does anyone know more about this?

Also, I would like to understand how to extend the two sets of frequencies to two sets of 8 frequencies, or 3 sets of 4 frequencies, such that it's possible to represent 64 keys. I'm not sure whether it would be easier to work with 3 sets of 4 versus 2 sets of 8 each. I guess the analysis would be easier if just using 2 sets of 8? There would be less intermodulation products to keep in mind …

Thanks!
B

EDIT: if it's too difficult to track back how this was achieved, i would like to know what the best way is to find 2 sets of frequencies, giving 64 unique combinations, and such that if you pick 2 frequencies from each of the 2 sets randomly, and mix the sine waves together, and put them through some non linear distortion, they are still easily detectable using goertzel, and the other frequencies are not "hit" (their bins remain at a normal background noise level)


texts: Why is the dB plot of a signal usually non-positive?

texts: Almost every decibel plot of a signal is usually below 0. Why is this the case? I recently just plotted a transfer function with a very high decibel plot above 0.


texts: A known signal, signal 1 got overlapped with an unknown signal. Likewise, signal 2 overlapped with another unknown signal. The problem I face now is how to classify the overlapped signals based on the signal source and label them as from signal 1 and signal 2 based on the available source training data of signal 1 and 2. Can you please help me in this?


texts:
In the A2DP specification, the SBC code uses a "polyphase filterbank" which has the formula
$$x[n] = X[m]\cos((m+0.5)(n+M/2)\pi/M).$$
I am confused to what this is? It seems to take in $N$ points and output $2N$ points in some sort of frequency domain. Is this a special type of DFT and can it be optimized with a radix-2 algorithm? I tried searching but could not find any websites that could help me.

Edit: Links - https://www.bluetooth.org/docman/handlers/downloaddoc.ashx?doc_id=8236 and bottom of Page 64


texts: How do I numerically calculate a function from its noisy gradient?

texts:
I have the model $\ s(x,y)=x^2+y^2, 0 \leq x \leq 1, 0 \leq y \leq 1 $.

Instead of observing the model directly I am observing the derivatives of the model + some noise (e):

$\ p(x,y)=s_x+e, q(x,y)=s_y+e$

From measurements of p(x,y and q(x,y) I want to estimate s(x). Say I know that s(0,0)=0.

According to the gradient theorem:
$\ s(x,y) = \int_{(0,0)}^{(x,y)} [s_x,s_y]d\vec{r}$

regardless along which path we integrate.

As a small experiment (in Matlab) I added normal distributed noise, N(0,1), to p=2x and q=2y. Then I integrated first along x followed by along y: SXY. Next I integrated first along y followed by along x: SYX.

The results show that the gradient theorem do not hold in this case (because of the noise):

The Root Mean Square Errors relative to the model are:

ErmsXY =
    0.1125
ErmsYX =
    0.0920

How can I find a better (less RMS error and smoother) estimate of s from p and q?

EDIT:

From what I read; using the curve integral is referred to as local integration.
There are also global integration methods where one instead tries to choose a S(x,y) that minimizes:

$\ \int_0^1\int_0^1 [|S_x - P|^2  + |S_y - Q|^2] dxdy$

Global integration methods are supposed to give better results when the gradient is noisy,
but how do I do this in practice?

EDIT 2:

One approach that I have used is this:

first we introduce linear derivation operators: $\ s_x = D_x*s,s_y = D_y*s$.

The result is the following linear equation system:

$\ D_x*s=p, D_y*s=q$

Next find a Least Square Error solution to these equations.
A LSE solution to these equations is supposed to be equivalent to minimizing the integral from above.
How can this be shown?

The results are good:

The RMS error is about 1/5 of that of SXY and SYX and the solution is also smoother.

However there are some drawbacks to this approach:

it is difficult to implement; must use central differences and "flatten" 2D s matrix into vector etc.
The derivation matrices are very large and sparse so they may consume a lot of RAM.

Another approach that seems potentially both simpler to code, less RAM consuming and faster is to use FFT. In the Fourier space these pdes become an algebraic equation.
This is known as the Frankot-Chellappa algorithm, but unfortunately I have not got it to work on my example data.


texts:
I'm trying to filter out bad face images (with artifacts) generated by StyleGAN.
Is there a good image quality measurement mechanism in python for doing this?
I've tried using BRISQUE but it is not showing good results (there are a lot of good quality images with higher scores then for bad quality images).

EDIT: These artifacts could be any kind of distortions, like eyes looking in a different directions, asymmetries (e.g. could be bad looking ear, assymetric eyes, ...), strange skin coloring, image messed up with the background like below image.

More images could be found here.


texts: Can error in digital signal be zero after filtering? I have looked up online to know about this but I could not find a solution.


texts: Difference between Uniform Linear Array (ULA) 3 dB beamwidth and bearing resolution

texts: Is there  a difference between linear array 3dB beam width and linear array bearing resolution ?


texts: Does Short-time Fourier transform impact quality of signal?

texts:
I am playing with spleeter to separate voice from an audio signal.
They go to the frequency domain by computing the Short-time Fourier transform to create a spectrogram. As there are some approximation (windowing and overlap of window) in this calculus, I was wondering if it hurts the quality of the signal when one computes stft and then inverse-stft.


texts:
I know that if a system is asymptotically stable then after infinite time its output  contains only characteristic of input signal and  it is also known as steady state and there are no transient component left in output, and if system is asymptotically unstable then whatever will be input is its output will become infinite after infinite time so we can say that it doesn't reach at steady state and it's transient components become infinitely large,
but what about transient and steady state in marginally stable systems?
Because in marginally stable systems even after infinite time it contains both inputs characteristics and systems characteristics so what exactly transient and steady state components will be?
I don't know whether I understand transient and steady state conditions properly
(whatever I explained above)  or it is just poorly defined in books?


texts: To implement successive interference cancellation (SIC) process, should the receiver know the power allocated to each symbol?

texts:
My question is as follows:

To implement successive interference cancellation (SIC) process in downlink non-orthogonal multiple access (NOMA), should receivers know all power allocated to each symbol?

Consider a single-cell network consisting of one base station (BS) and $N$ users.
Let $x_i$ be the symbol for user $i$, $h_i$ be the complex-valued channel gain between user $i$ and the BS, $p_i$ be the power allocated to symbol $x_i$.
Then, the transmitted signal is given by
$$
x = \sum_{i=1}^N \sqrt{p_i} x_i.
$$
The corresponding received signal at user $k$ is given by
$$
y = h_k \sum_{i=1}^N \sqrt{p_i} x_i + n_k,
$$
where $n_k\sim\mathcal{CN}(0,\sigma^2)$ is the additive white Gaussian noise of user $k$.
Assuming that $\lvert h_1 \rvert^2 \le \lvert h_2 \rvert^2 \le \cdots \le \lvert h_N \rvert^2$, user $k$ decodes its own symbol $x_k$ according to the following process:

Decode $x_1$ by treating $\sum_{i=2}^N h_k\sqrt{p_i}x_i + n_k$ as noise.
Subtract $\hat{x}_1 \triangleq h_k\sqrt{p_1}x_1$ from $y$.
Decode $x_2$ by treating $\sum_{i=3}^N h_k\sqrt{p_i}x_i + n_k$ as noise.
Subtract $\hat{x}_2 \triangleq h_k\sqrt{p_2}x_2$ from $y-\hat{x}_1$.
Repeat to (i) decode $x_l$ by treating $\sum_{i=l+1}^N h_k\sqrt{p_i}x_i+n_k$ as noise, and (ii) subtract $\hat{x_l} \triangleq h_k\sqrt{p_l}x_l$ from $y-\sum_{i=1}^{l-1}\hat{x}_i$ until $x_k$ is decoded.

In my opinion, to decode its own symbol $x_k$, user $k$ should know $p_l$ for any other user $l$ such that $\lvert h_l \rvert^2 < \lvert h_k \rvert^2$.
Is this opinion true? If this is true, the BS should inform the user about the power allocated to each symbol. Nevertheless, is NOMA useful?


texts: Question regarding match filter implementation in receiver

texts:
For example, if you have a QAM 16, which mean there are 16 symbols.  Each symbol will have different level of the QAM constellation.  Each symbol will have a corresponding match filter.  At the receiver, there is a bank of match filter to match with each received symbol.
The receiver will use threshold to detect which has the highest threshold of that matched filter array.  My question is, let's say you have two symbols: one with high level and one with lower level.  Is it guaranteed that only the symbol that got "matched" will have the highest threshold level after the matched filter output?
I mean, can the high level symbol have a higher output after the matched filter even if it is not matched vs. a lower level symbol that is matched?


texts: Why do people sometimes write $H(e^{j\omega})$, some others use $H(\omega)$ and even others use $H(f)$ with $\omega=2\pi f$ to describe the frequency response of a filter or a spectrum? Is there any "real" difference in meaning or is it just convention?


texts:
I've been learning to implement an IIR filter for a synthesizer project. It is now more or less stable when the parameters are constant, but ultimately I want the parameters to be arbitrarily modulated.

I figured it would be too slow to compute new recursion coefficients every sample, so instead I'm computing them every "chunk" (a piece of 256 or so samples which the synthesizer produces in a batch and then sends to the audio driver). This brings me to my question:

How should I approach hotswapping the recursion coefficients? Or is there a more correct way of modulating IIR filters?

I am curently linearly interpolating the output of the model generated in the last chunk with the one generated in the current. (at the beginning of the chunk, the output is 100% from the old one, and at the end 100% from the new one) This isn't enough, since the new filter needs to be initialized somehow. I have tried some things with various degrees of success:

The new filter is initialized with signal history pulled from the old one.

This sounds the best, but when too heavy modulation is applied (for example, the cutoff frequency changes too much at once), the new filter might initially output an uncontrollable spike which then peaks. This mostly happens along with a discontinuity in the input signal.

Here's a 60 Hz saw wave LP-filtered using a cutoff oscillating at 20 Hz. (the filter is configured to have two biquads). Such a spike never occurs when the cutoff frequency is kept constant.

The filter history is initialized in a stable state. (e.g. constant signal in and out at the intensity of the last seen sample for LP, or plain zero)

This causes some audible rattle when modulating, no matter how I try to fade in the new filter gradually. The new filter never quite reaches the old one's state.

All ideas and advice are welcome.

Edit:

I have implemented the biquads in a lattice topology, as per Robert's suggestion. Modulation now works great, but raising the cutoff frequency too high (or too low for a high pass filter) breaks something. This should be a single biquad low pass filter with cutoff $\frac{3}{4}\pi$:

The direct form parameters, from which the lattice is derived, work just fine. It is certainly possible, that there is some error in my equations for the conversion. I also noticed that the shape of the frequency response isn't quite as steep as with the direct form. I am suspecting it has something to do with the "tap" coefficients..

What could be going on now? Perhaps I should make a new question with details about my implementation, as the original issue is now more or less solved.


texts: What is the right algorithm to detect segmentations of a line chart?

texts:
To be concrete, given 2D numerical data as is shown as line plots below. There are peaks on a background average movement (with small vibrations). We want to find the values of pairs (x1, x2) if those peaks drops down to average; or (x1) only if the line doesn't back to the average.

There are thousands of such 2D data.

What is the right statistic or machine learning algorithm to find x1 and x2 above without plotting?


texts: Are all images 2d or also we have 3d images?

texts: Are all the images that we encounter/use in image processing are 2D or do we also have 3D images? If also 3D, then Video will be treated as a 2D image or 3D image?


texts:
I am trying to design a low pass filter to apply before resampling the sensor's output at a lower rate, however I am not quite sure how to choose between the IIR and FIR fitlers. Matlab by default, filters the data with an 8th order Chebyshev Type I lowpass filter with
cutoff frequency .8*(Fs/2)/R, before resampling. but how can I choose the best low pass filter type ?


texts:
I have a rectangular matrix S which is the amplitude of S-transform and a matrix G whose arrays are zero and one.
G and S have the same size. Each row corresponds to time ti and each column to frequency fj.
Now, I have to calculate the convolution of G*S in time frequency domain such the result should be the same size of S.
How can I do that?
Do you know any MATLAB function which does this calculation?
Thank you so much


texts: If a time-series has odd number of samples does it have no energy at Nyquist frequency?

texts:
Suppose I have real time series A with n samples and time-spacing dt and I want to analyze its frequency content.
Af = fft(A)

If dt=1 Nyquist frequency is 0.5.
According to both Numpy and Matlab manuals, the frequency vector is defined as
f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (dt*n)   if n is even

f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (dt*n)   if n is odd

If n=9 (odd), f= [-0.44444 -0.33333 -0.22222 -0.11111  0.  0.11111 0.22222  0.33333  0.44444]
If n=8 (even), f= [-0.5   -0.375 -0.25  -0.125  0.     0.125  0.25   0.375]

Note the maximum frequency for odd number of samples is 0.44444 Hz not 0.5 Hz. Does this mean a time series with odd number of samples has no energy at Nyquist frequency ?


texts:
This is a Question asked in  IISC ( Indian Institute of Science,Bangalore,India) interview for MS admission.

What is the $\mathcal{Z}$-transform of  $\dfrac 1{n^2}$ ?


texts:
Why does cel-shading sometimes look similar to rotoscoping?

Rotoscoping: A Scanner Darkly:

Cel-shading: Wikipedia: Cel-shading, Teapot Example:


texts: Properties of Spectral Transformations - Allocation (decomposition into even and odd part)

texts:
I am trying to understand the Allocation property of Spectral Transformations. I can't.

I know that every function can be separated into an even part and into an odd part.

My problem is understanding proof for that:

$$s(t)=s_e(t)+s_o(t)$$
where:

$s_e(t)$ - even part of function
$s_o(t)$ - odd part of function

The problem is that in next step at every site I looked so far, it goes something like this:

$$s_e(t) =\dfrac{s(t)+s(-t)}{2}$$

Why did they do this? Where did $\frac{1}{2}$ come from, and why is $s(t)+s(-t)$

Can you please explain it to me, or point me to a good resource where I could read more about it?

Thanks!!!


texts: How can I increase image size by zero padding?

texts:
For increasing the image size from 8X8 to 32X32, I tried to do zero padding. The simplest concept of zero padding in Java for me seems to be the code as shown below in Java such that I can shift each point by 4 points while rescaling.

for (int i = 0; i < imageWidth; i += 4)
        for (int j = 0; j < imageHeight; j += 4)
            pspectrum[i][j] = powerspectrum[(i + 1)/4][(j + 1)/4];

Now, what I am mainly concerned is that is there any more accurate method and how this helps with my calculation.


texts: Find points that are mapped to the original grid with image rotation by a non multiple of 90 degrees

texts:
I would like to find the point that are "not interpolated" when I rotate an image with e.g. Matlab's imrotate() by an angle which require interpolation e.g. 6 degrees, i.e. not a multiple of 90 degrees.
In other words, I would like to perform rotation by an arbitrary angle but WITHOUT any interpolation. (i.e., I expect to find "holes'' where the fractional coordinates are computed by the multiplication with the rotation matrix.)
The idea I had was to use the points which would be mapped to an integer coordinate after rotation by the rotation matrix [cos(t) -sin(t);sin(t) cos(t)] but that does not seem to work since there is always some approximation error in a numerical analysis oriented language like Matlab, but in computers in general of course.
Any ideas/functions?


texts: Introduction to speech processing with code and examples?

texts:
Can anyone suggest an introduction to speech processing,
class notes or books
accompanied by programs that can be run to see what the algorithms do ?
For example, in Burrus et al., Computer-Based Exercises for Signal Processing Using Matlab,
20 of the 400 pages are on speech processing. But it's from 1994;
there must be other hands-on books ?
Books with reviews like the following may have their uses, but are not what I'm looking for:

  I found this book comprehensive but incomprehensible, primarily because of the
  lack of real-world examples. It would benefit from a series of programming exercises
  with training sets

(This question may be a tall order — not many people ar good at theory and practice and explaining.)

(Should this be a community wiki ?)


texts:
Ok, so I'm trying to compare two different speech signals and I have come into a problem. Here goes:

I have split the signal into blocks, and I have computed the MFCC coefficients of each block. I then use a DTW algorithm to compare the (inputted) signal to the training signal.

Here is the data:

The training data is the MFCC values of a signal of someone saying "No". If I have the inputted signal of someone saying no the values then give:

  13.9462,
  55.8784,
  38.3383,
  29.9468,
  32.7136,
  24.8893,
  34.0734,
  24.3645,
  39.329,
  20.4847,
  31.1939,
  29.8841,
  25.7655,
  28.222,
  23.1643,
  33.4366,
  .....,
  ......

Where the max data for this dataset is: 97.4834

Here is the data for someone saying a different word:

  28.6696,
  65.8777,
  44.2725,
  31.6083,
  42.6541,
  38.4104,
  26.6311,
  34.9188,
  37.2065,
  25.2479,
  41.5969,
  54.2681,
  37.0685,
  26.2073,
  33.9836,
  38.7847,
  28.3622,
  67.8788,
  74.9075,
  .....,
  .....
  Where the max value in this dataset is: 97.5609

My guess (reading through the algorithm) is that the smallest difference should be the correct result. However, I do not know if it is possible to either:

1) Calculate the MAX value for each dataset and then see which value is smaller

2) Get the first variable and compare it.

Does anyone have any suggestions to how I would compare these kind of values?

Thank you


texts: Output of a linear time-invariant(LTI) system

texts: I am very confused about one of the questions I received during an exam. How do I solve this question?


texts: What's the relationship between the parameters of this noise generator and the gaussian curve it produces in the frequency domain?

texts:
This algorithm produces a sequence $y$ having a magnitude spectrum that resembles a normal distribution, peaking at frequency $\omega$:

$$b[n]=(1-\sigma)b[n-1]+\frac{\sigma}{\omega} C$$
$$p[n]=p[n-1]+\omega(1+b[n])$$
$$y[n]=cos(p)$$

where $b[n]$ is a leaky integrator having feedback $1-\sigma$ and taking white noise $C$ as input, $p[n]$ is a phase accumulator having frequency $\omega$ modulated by brownian noise $b[n]$, and y[n] is the final signal, cosine at phase $p[n]$.

Python function:

def signal(f, w, n):
    """
    Arguments
    f - peak frequency in radians/sample
    w - bandwidth of frequency content
    n - length
    Returns
    sig - list of samples

    """
    br = 0
    p = 0
    sig = [0] * n
    for i in range(n):
        # br - leaky integrator having rate proportional to bandwidth,
        # integrating white noise at magnitude inversely proportional
        # to frequency.

        br = (1 - w ) * br + w * random.uniform(-1, 1) / max([0.001, f])

        # p - phase increments at frequency f modulated by brownian noise.
        p = p + (1 + br) * f

        sig[i] = math.cos(p)

    return sig

To illustrate:

averaged magnitude spectrum at swept frequency from [0,π)

averaged magnitude spectrum with swept bandwidth from [0,0.5)

fourier transform of a sample sequence, hue indicates phase

wav

audio sample at 1760hz with $\sigma=0.005$

wav

averaged eight times

There are undesirable characteristics:

The amplitude of the peak frequency has an unknown relationship to $\sigma$, the second animation shows the sharp curved response
$\omega$ is not a perfect predictor of peak frequency, $\sigma$
allows noise to perturb the phase towards a range of
$[0\omega,2\omega]$ from FM

Short of using optimisation and curve-fitting, how can the frequency content be predicted (e.g. by the gaussian function) and controlled?

Any additional insight appreciated


texts:
Here is a 10 seconds-long 440hz sine wave
 normalized at $0\textrm{ dBFS}$.

When computing the STFT (with the code below) of this audio file, I noticed that max(abs(STFT)) is around 248.33.  (more generally, it seems to be approximately fftsize/4 for this particular file).

Even after doing magdB = 20 * math.log10(abs(STFT)) to have decibel values, we get a max of 47.9 dB.

But this $47.9\textrm{ dB}$ doesn't mean anything! We would like to have ~ $0\textrm{ dB}$ instead because this audio only contains 1 sine wave component (1 "harmonic") at $0\textrm{ dB}$ volume.

Audio tools display it properly, here is the spectrum display of this sine wave:

Question: how to have an absolute, canonical $\rm dB$ values from a FFT, that makes that a pure $0\textrm{ dBFS}$ sinewave has a peak of $0\textrm{ dB}$ in the spectrum display ?
More generally, is there a canonical way to go from FFT values to $\textrm{dB}$ in order to display a spectrum analysis?

NOTE: I know that $\textrm{dB}$ are for ratio between 2 things, so all I mentioned here is referred to $\textrm{dBFS}$ (full-scale)

import scipy, numpy as np
import scipy.io.wavfile as wavfile

def stft(x, fftsize=1024, overlap=4):
    hop = fftsize / overlap
    w = scipy.hanning(fftsize)
    return np.array([np.fft.rfft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])


texts: Does 46 dB gain of a filter for a frequency imply 200 times more amplitude?

texts: Consider we have a IIR filter or order 2 and it has gain of say 46.02 dB for frequency $f_0$. Now if we apply a time domain signal of frequency $f_0$ to the filter for a long time, can we expect to get 200 times higher amplitude in the output of the filter comapre to the input of the filter? The 46.02 db and 200 comes from the fact that $20\log_{10}(200)= 46.0206$.


texts: What are the differences between the various FFT algorithms?

texts:
Wikipedia lists a number of different FFT algorithms:

Cooley–Tukey FFT algorithm, Prime-factor FFT algorithm, Bruun's FFT algorithm, Rader's FFT algorithm, and Bluestein's FFT algorithm

What are the pros and cons of each algorithms? (think about: computational expense, parallelizability, ease of implementation, etc.)
And which ones are preferred for which applications?


texts: How to set sampling frequency of an acceleration data set

texts:
I have a set of accelerometer readings in $X, Y, Z$ axes obtained from an android based smart-phone. The data in $(X, Y, Z)$ was recorded at different time stamps and there is no uniform time period for recording the data.

My data set looks as  timestamp $X, Y, Z$. First, I did a filtering (using a low-pass) on this data and now want to perform $\textrm{FFT}$ on a timewindow of 1 min or may be a window containing 250 samples (not very sure about the window length). I am using this FastFourierTransform class of apache commons in Java.

I am getting $\textrm{FFT}$ magnitude but I was wondering how do I know the corresponding frequency of each $\textrm{FFT}$ mag. I know corresponding frequency would be $(n\cdot F_s)/N$ where $n$ is the bin number or data point index, $F_s$ is sampling frequency and $N$ is number of input data points over the window. Now my question is how do I know $F_s$ for a given set of data say for example input data in an array as $\begin{bmatrix}1&2&3&4&5&6& 7\end{bmatrix}$ over a window of size $N_t$ millisec where $N_t=(\textrm{lastTimeStamp}-\textrm{firstTimeStamp})$

Is it $\displaystyle \frac{7}{(N_t\cdot 10^{-2})} \textrm{Hz}$?

This also immediately follows by doing so if this agrees with nyquist sampling frequency? and this says $F_s\ge2F_{max}$ in the signal and I don't know the maximum frequency over the time window $N_t$.


texts: How to do a Fourier transform on input from a microphone and obtain the average amplitude of a certain frequency range

texts:
I'm building a sound visualizer with LED strips controlled by python, and I'm trying to learn how to use pyAudio to perform rapid Fourier transforms on live mic input and obtain the average amplitude of a certain range of frequencies.

How do I open a sound stream and consistently update the current amplitude (say, a variable called "amplitude") of a given frequency range with the live FFT data? I have set up all the hardware and all the basic pyAudio functions work so far.


texts: Why doesn't a sudden loud noise sound high pitched?

texts:
The jumpiness or high change of a signal is due to the higher frequency components in the signal. So if I have a sound signal that increases suddenly in amplitude, why doesn't the signal sound very high pitched at that point?

For example, if I suddenly knock on a door or hit something, the pitch doesn't sound very high, but the magnitude of the sound increases from 0 to a number very quickly.


texts:
I would like to implement a HOG descriptor in C++.
I found an implementation of this code here.

An extract of the code follows:

//extract HOG features from a group of key points
void HOGExtractor::extract(ublas::matrix<float> &dscr_list, const
ublas::matrix<int> &bbox_list, const IntegralHistogram &inthist)
{
    unsigned int length = _param.xgrid * _param.ygrid * inthist.dirnum();
        if (dscr_list.size1() != bbox_list.size1() || dscr_list.size2() !=
length)
    throw std::runtime_error( boost::str(
                boost::format("Allocate proper memory for dscr_list [%1% %2%] != [%3% %4%]") % dscr_list.size1() % dscr_list.size2() % dscr_list.size2() % length) );

ublas::matrix<float>::iterator1 idscr_list;
ublas::matrix<int>::const_iterator1 ibbox_list;
ublas::matrix<int>::const_iterator1 ibbox_list_end(bbox_list.end1());
ublas::vector<int> bbox_grids( _param.xgrid * _param.ygrid * 4 );
for ( ibbox_list = bbox_list.begin1(), idscr_list = dscr_list.begin1(); ibbox_list != ibbox_list_end; ++ibbox_list, ++idscr_list)
{
    ublas::matrix<int>::const_iterator2 ibbox = ibbox_list.begin();
    float x0 = std::max( *ibbox++, 0 );
    float y0 = std::max( *ibbox++, 0 );
    float x1 = std::min( *ibbox++, inthist.width() - 1 );
    float y1 = std::min( *ibbox++, inthist.height() - 1 );
    if ( x1 - x0 < _param.xgrid || y1 - y0 < _param.ygrid )
        throw std::runtime_error( boost::str(boost::format("Grid number should be less than bbox width/height [%1% %2% %3% %4%]/[%5% %6%]") % x0 % y0 % x1 % y1 % _param.xgrid % _param.ygrid) );

    float xstep = (x1 - x0 + 1) / _param.xgrid;
    float ystep = (y1 - y0 + 1) / _param.ygrid;
    for (int i = 0; i < _param.xgrid; ++i)
        for (int j = 0; j < _param.ygrid; ++j)
        {
            bbox_grids( (i * _param.ygrid + j)*4 + 0) = static_cast<int>(x0 + xstep*i);    //x0
            bbox_grids( (i * _param.ygrid + j)*4 + 1) = static_cast<int>(y0 + ystep*j);    //y0
            bbox_grids( (i * _param.ygrid + j)*4 + 2) = static_cast<int>(x0 + xstep*(i + 1) - 1);    //x1
            bbox_grids( (i * _param.ygrid + j)*4 + 3) = static_cast<int>(y0 + ystep*(j + 1) - 1);    //y1
        }
    inthist.get_hist< ublas::matrix<float>::iterator2, ublas::vector<int>::const_iterator>( idscr_list.begin(), bbox_grids.begin(), bbox_grids.end(), _param.normalize);
}

}

The "key points" it is referred to here are (if I'm not mistaken) the coordinates of the box containing an object we want to extract the features.

However, I am wondering how I can get these coordinates in a image? (Given an image and an object I want to extract the features, how can I know the values of x0, x1, y0, y1?)


texts: Non-coherent receiver for DBPSK: Implementation details

texts:
I am implementing a non-coherent receiver for DBPSK.
I have a computer science background with little knowledge of DSP so my questions here may seem obvious.
The design of the receiver I am trying to implement is the following:
Design of DPSK receiver

And what I am doing right now is:

On the tx side

Differentially encoding my bit sequence
BPSK modulating it
Applying a high-pass filter

On the rx side

Applying a band-pass filter around the carrier frequency
Delaying (by one bit worth of data) and multiplying the data by itself
Applying a low-pass filter
Integrating over one bit worth of data (i.e. sum the sample values)
Making the decision depending on the sign of my integral

By 'data' I mean samples at 44.1kHz.
When I simulate this on Octave, it works (i.e. I can detect the bit sequence even with high AWGN).
The problem is I am trying to move to an infinitely long data sequence (i.e. I do not know the exact moment where the tx began).
So here are my questions:

How can I achieve time and frame synchronization ? (I tried with a barker code but I am not sure I am getting it right)
This is the thing that's most confusing to me. My guess is you'd have to get time synced before the integration block, and frame synced after the decision block.
How can I express the bandwidth in function of carrier frequency f and bit rate ?
Does the relation between bit rate and carrier frequency have any influence on the receiver ?
What are the optimal values for the filters ? (number of taps)
How do I calculate Eb/N0 for in my case ? (I want to compare my BER curve with others to see how good is my receiver)


texts: Determining type of filter given its equation

texts:
Given a general filter equation, how can one tell the type of filter that the same equation represents? Meaning how can I know if the filter is Low/High/Band Pass etc...?

For example, the following equation:

$ y_{n} = x_{n} + x_{n-2} $

Represents a "Band Stop" filter, but why?


texts:
I'm trying to understand how diode circuits are implemented in wave digital filters, particularly for clippers. The research papers and other sources I've looked at use the equation
$$I(V) = 2I_s \sinh\left(\frac{V}{V_t}\right)$$
for two reverse-polarity diodes in parallel, then transform to the wave domain using
\begin{align}
V &= \frac{a+b}{2}\\
I &= \frac{a-b}{2R_p}
\end{align}
Where $R_p$ is the port resistance, $a$ is the incident wave and $b$ is the reflected wave:
$$
\frac{a-b}{2R_p} = 2I_s\sinh\left(\frac{a+b}{2V_t}\right)$$
I'm not the best with WDFs, but my understanding is that the incident wave a is propagated from the rest of the circuit up to the root (the nonlinearity, in this case the diodes), and then b is computed and reflected back down through the circuit. So a is the given in the above equation. But this is still not a trivial equation to solve, and has multiple solutions if I'm not mistaken.
None of the resources I've found actually show the solution to this equation, and instead just express it as b = f(a). This paper briefly mentions using a Lambert function?
This paper is the closest I've found to being helpful, but even here the conditions seem rather arbitrary? (What are G0, Gv or Gn for a diode??)
The last source I've looked at is this one which has some MATLAB code for a circuit with a singular diode. Don't know if anybody knows more about this stuff than I do, but if you have any pointers it would help a lot.
EDIT: after re-visiting the MATLAB code it looks like they're using a combination of the relationship between b and a for linear resistances:
b = a(R-Rp)/(R+Rp)
and substituting R for the nonlinear diode resistance, which according to this source is
R = Vt/I
I'm not 100% sure that's what's happening there, but it's the closest I've found to a tangible solution thus far


texts: What are the applications of Fourier Transform Power Spectrum?

texts:
What are the main applications of Fourier Transform Power Spectrum? Can somebody clarify with some examples.

Also, what are the advantages of using Fourier Transform Power Spectrum over Fourier Transform?

What is the advantage of Fourier Transform Power Spectrum in case of Image Processing?

Thank you so much.


texts: How to compare the similarity of 2 transfer functions

texts:
I'm working on a power systems where I have 2 linearized transfer functions. The first transfer function corresponds to the no-load case and the second transfer function corresponds to the full-load case.

I was wondering if there are metrics to compare 2 transfer functions, for example how similar are their respective poles and zeroes.

I can play on some factors, like inductance and switching frequency and I get the feeling that increasing the inductance value too much makes the 2 transfer functions more dissimilar.

Regards


texts:
The idea of OFDM is to spread data across multiple overlapping subcarriers. As they are orthogonal it is still possible to demodulate them. That means, all other subcarriers are 0 at the evaluation point of a specific subcarrier.

Orthogonality is achieved by generating a rectangular subcarrier pulse. By applying a IFFT a $\sin(x)/x$ emerges which allows orthogonality.

Another requirement for orthogonality concerns the subcarrier spacing. It must be equal to the reciprocal of the symbol period. Howver, I do not understand why. I thought it is most important to create these $\sin(x)/x$ that are shifted across the x axis on an unique position for each subcarrier.

Question:
Why would orthogonality break if the subcarrier spacing was not reciprocal of the symbol period?


texts:
An FIR linear phase filter has unit sample response $h[n]$ that is real with $h[n]=0$ for $n < 0$ and $n > 7$.  If $h[0]=1$ and the system function has a zero at $z=0.4e^{j\pi/3}$ and a zero at $z=3$, what is $H(z)$?

Solution says, that it has $7$ zeros and here's the system equation:

$$
H(z)=A(1-0.4z^{-1}+0.16z^{-2})(0.16-0.4z^{-1}+z^{-2})\cdot \\(1-3z^{-1})(1-\frac{1}{3}z^{-1} )
$$

My question is:  How does this system function $H(z)$ have $7$ zeros?

When I count them I see:

2 zeros for: $(1-0.4z^{-1}+0.16z^{-2})$

2 zeros for: $(0.16-0.4z^{-1}+z^{-2})$

2 zeros for: $(1-3z^{-1})(1-\frac{1}{3}z^{-1} )$

6 zeros total.

The book says it has 7 zeros. however, i'm not seeing the 7th zero. How to find it?


texts: Does "improper" imply that a system cannot be stable and causal?

texts:
This answer and the comments in it made me wonder whether the following statement is true:

  If a transfer function is improper, then that system cannot be causal and stable at the same time.

I had thought that this was true for a while. But the other day I wondered why. For example, the transfer function

$$H(s)=\frac{s^2}{s+1}$$

is improper, but the ROC $\{s:\mathrm{Re}(s)>-1\}$ would make it stable (it contains the imaginary axis) and causal (it consists of a left-sided plane).

So how are causality and stability related in improper systems?


texts:
I'm trying to process an 1D signal containing a waveform of 100000 pixel.
As it was acquired with a laser scanning system, and therefore it is corrupted by some annoying speckle noise.
I tried to prefilter it with a well-known type of non-linear image processing filter such as Lee and Kuan filters, obviously modified for the 1D nature of the signal.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2784081/ (see citation 5,7)

However, the resulting signal is identical to the result of a "naive" mean filter with the same support size (e.g. 15 px).

Is it an "expected" behaviour due to the 1D processing or am I probably writing wrong code?
Can you link me some working code for these algorithms or some other speckle reducing ones?
I found some source code on codeforge, but I think the Lee and Kuan algorithm are not correctly implemented.

Thank you very much in advance.

EDIT:
added some images of the signal:

1) original

2) mean 1x15

3) mean 1x15 with superimposed kuan filter 1x15 with input parameters  nm = 1.0 and nv = 0.1

4) mean 1x15 with superimposed kuan filter 1x15 with input parameters  nm = 1.0 and nv = 0.2

5) mean 1x15 with superimposed kuan filter 1x15 with input parameters  nm = 1.0 and nv = 0.3

We can see in 3,4,5 that there're no differences by mean and kuan: the mean signal is exactly coinciding with the kuan filter one.

Note that at the pixel 50569, where there is an exampe of the speckle I want to reduce, the local mean is 121.666 and local variance is 16.755. Therefore it can be found that an estimation of the noise variance is 0.13771, then I tried the parameters according to said observation.
Further, note that at pixel 50652 the signal is lower and therefore de speckle effect is lower too. There, the local mean is 40.933 and local variance is 4.195 and the ratio is 0.102, therefore in the order of the ratio on an high peak.

Thank you again for your attention


texts: Under what conditions is there a one-to-one mapping between continuous-time and discrete-time signals?

texts:
As the sampling theorem dictates that the uniform sampling frequency must be at least twice the maximum frequency present in the bandlimited signal (Nyquist rate), a question arises about the uniqueness of this representation of the analog signal to digital one.
Apart from practical errors related to the truncation of the infinitely-long energy signal when discretized to digital, can we find a UNIQUE mapping between the two signal forms? If yes, then what are the conditions that limit or control this unique representation?
Note: I am not assuming by default that the sampling topology should be uniform. Random sampling also could be considered as well.


texts: What is the name of this digital low pass filter?

texts:
I found this digital filter in code I am working on. It is a low pass filter. In the code it is called an "alpha filter", but it is not the same as the alpha filter mentioned here.

I post the relevant code below:

float alpha = 0.7;
float prev = 0.0;

float filter(float sample) {
    float filtered = (sample * alpha) + (prev * (1.0 - alpha));
    prev = sample;
    return filtered;
}

It looks like alpha should be in range $[0, 1]$. By increasing alpha, the filter tracks the measurement with less delay, but lets more noise through.

I am tasked with writing library functions, and I would like to give this thing its proper name, and hopefully link to a Wikipedia article in the doc tags.


texts: Memory requirements for FFT on STM32F103C8

texts:
I have a limited system and would like to implement FFT within STM32F103C8 without any extra memory buffers.

So I want to know how many memories are needed to implement if I have 2592x1944x8bit size's one image?

Actually, I want to have a process such as

Origial image ---> FFT ---> Blur ---> IFFT ---> Modified image
The blur(Gaussian) size is quite a bit huge size. 60x60~90x90 variable.

What is the memory requirements for FFT on STM32F103C8 ?

The reason of why I pick this way up is 60x60~90x90 gaussian blur need a huge memory to process to blur. and  I found FFT blur processing from google, so this is the reason.


texts: Extract breaths from respiratory signal, i.e. peak detection

texts:
I am trying to extract breaths from a respiratory signal, which is essentially peak detection. See the example below for what such a signal looks like. This particular signal looks rather clean, however, lower amplitude noise is very common (both high and low frequencies).

I have tried the ones mentioned here, however, those seem to require that the height of the peak can be described by its y-value. As you can see below, the signal is very prone to baseline wander, and the top of a peak can very well fall below 0.

Any suggestions for a peak detection algorithm?


texts:
I have one image. I want to segment it by simple approach. Could you suggest to me the best method to segment it? I think one of them is gradient or FCM. But I am not sure. Please see  my given image file here
and segmented image file below


texts: pole/zero locations for real and imaginary signal

texts:
In Z-Transform,

For a real signal,
 $x(n)$ =$x^*(n)$ .
Taking Z-transform on both sides,

$X(z)$=$X^*(z^*)$ ,  which gives certain pole/zero condition

similarly for a purely imaginary signal            ,

$x(n)$ =$-x^*(n)$ , if we take Z-Transform on both sides , $X(z)$=$-X^*(z^*)$

which gives similar pole/zero condition as that of real signal(poles and zeros are not effected by negative sign).

But  My professor says if the poles/zeros are symmetric with respect to real axis, then it is a real signal. But i am getting same condition for both real and imaginary signal.

Can anyone help me? (examples might be helpful)


texts:
I am reading a paper on SIFT by Lowe.

In it, gradient magnitudes and orientations are calculated using these formulae:

$$$$
I don't know why m(x,y) is calculated above. I heard that gradient magnitude is usually calculated like the following.

$$\sqrt{f_x ^2 + f_y^2}
$$


texts: calculating DC value from positive frequency with its DC undefined

texts: Is there a way of calculating DC value from just the positive frequency that has its DC undefined ? Simply if we remove the DC and the negative frequency components of the signal is there a way to reconstruct it from just the positive frequency values.


texts: Programming Wavelets for Audio Identification

texts:
How exactly is a wavelet used digitally?

Wikipedia states:

  For example, a wavelet could be created to have a frequency of Middle
  C and a short duration of roughly a 32nd note. If this wavelet was to
  be convolved with a signal created from the recording of a song, then
  the resulting signal would be useful for determining when the Middle C
  note was being played in the song.

Would this be a data structure holding e.g {sampleNumber, frequency} pairs?

If a wavelet is an array of these pairs, how is it applied to the audio data?

How does this wavelet apply to the analysis when using an FFT?

What is actually being compared to identify the signal?


texts:
This question stems from an article,

"Entropy estimation of symbol sequences"  download link

where the abstract mentions the need for using symbols in information theory.

  Random chains of symbols (zeros and ones or characters) $s_1, s_2,\ldots,$ drawn
  from some finite alphabet
  appear in practically all
  sciences. Examples include spins in one-dimensional mag-
  nets, written texts, DNA sequences (here symbols are alphabets), geological records of
  the orientation of the magnetic field of the earth, and bits
  in the storage and transmission of digital data.

The paper presents technique for finding the entropy of symbols. My Question is that if entropy of symbols is important then what about the entropy of real valued numbers?
From the perspective of entropy, is there a proof or some information in text where it provides information that the entropy of random symbols is greater or lesser than the  entropy of real numbers.

Q1: In other words, what would be the entropy of random variables having binary values and entropy of real valued random variables?  A fair coin takes 2 values -0/1 with equal probability and its entropy is 1 bit. So, is 1 bit the maximum entropy? Would entropy of binary random variables be less or more than entropy of real valued random variables and how to show?

So, I would like to know if there is a proof or material from which it can be inferred if the entropy of characters is lesser/ greater or same with  the entropy of real numbers.

Q2: I want to know why symbols are used in coding and other applications such as storage from the viewpoint of entropy and information theory.


texts: Recovering time domain signal from sum of independent and identically distributed signals

texts:
Background: I observe a sample of a variable z that is the sum of two independent and identically distributed variables x and y. I'm trying to recover the distribution of x, y (call it f) from the distribution of z (call it g), under the assumption that f is symmetric about zero. According to Horowitz and Markatou (1996)  we have that the Fourier Transform of f is equal to sqrt(|G|), where G is the Fourier transform of g.
Example:
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde, laplace

# sample size
size = 10000
# Number of points to preform FFT on
N = 501
# Scale of the laplace rvs
scale = 3.0

# Test deconvolution
laplace_f = laplace(scale=scale)
x = laplace_f.rvs(size=size)
y = laplace_f.rvs(size=size)
z = x + y
t = np.linspace(-4 * scale, 4 * scale, size)
laplace_pdf = laplace_f.pdf(t)
t2 = np.linspace(-4 * scale, 4 * scale, N)

# Get density from z. Kind of cheating using gaussian
z_density = gaussian_kde(z)(t2)
z_density = (z_density + z_density[::-1]) / 2
z_density_half = z_density[:((N - 1) // 2) + 1]
ft_z_density = np.fft.hfft(z_density_half)
inv_fz_density = np.fft.ihfft(np.sqrt(np.abs(ft_z_density)))
inv_fz_density = np.r_[inv_fz_density, inv_fz_density[::-1][:-1]]
f_deconv_shifted = np.real(np.fft.fftshift(inv_fz_density))
f_deconv = np.real(inv_fz_density)

# Normalize to be a pdf
f_deconv_shifted /= f_deconv_shifted.mean()
f_deconv /= f_deconv.mean()

# Plot
plt.subplot(221)
plt.plot(t, laplace_pdf)
plt.title('laplace pdf')

plt.subplot(222)
plt.plot(t2, z_density)
plt.title("z density")

plt.subplot(223)
plt.plot(t2, f_deconv_shifted)
plt.title('Deconvolved with shift')

plt.subplot(224)
plt.plot(t2, f_deconv)
plt.title('Deconvolved without shift')

plt.tight_layout()
plt.show()

Which results in

Issue: there's clearly something wrong here. I don't think I should need the shift, yet the shifted pdf seems to be closer to the truth. I suspect it has something to do with the domain of the IFFT changing with the sqrt(abs()) operation, but I'm really not sure.


texts:
I have a serious doubt about Autocorrelation and AutoConvolution.
My understanding is,

Autocorrelation

Autocorrelation is nothing but correlation of the signal and it's delayed version. Example : Autocorrelation of 2 functions f(t) & f(t+1).

Autoconvolution

Autoconvolution is combination of  the two signals where the output signal will be a modified version of a signal in term of another signal with every point/unit.

My Question

Are two signals are the same if their auto-correlation functions are the same?
Are two signals are the same if their auto-convolution functions are the same?


texts:
It is easy to define an (ideal) LTI system that would have an infinite number of poles - for instance, if the transfer function is
$$
H(z)=\frac{1}{\cos(z)-1}
$$
However, this would only define a countably infinite set of poles.
I am curious: are there systems with uncountably infinite sets of poles?
I would guess that the very definition of a pole would prohibit a function from having a contiguous region of poles, but what about a hypothetical transfer function that would have its poles defined on something like Cantor set?


texts: How do you calculate the group delay in samples of a Biquad Filter at a given frequency from the coefficients?

texts:
If I have a BiQuad all-pass filter with a given set of coefficients such as:

g: 0.903228
a1: -1.88142
a2: 0.903228
B1: -2.08299
B2: 1.10714

And my operation per sample to process the filter is:

    inputGained_2 = inputGained_1;
    inputGained_1 = inputGained;
    inputGained = input * g;
    output_2 = output_1;
    output_1 = output;

    output = inputGained + (B1 * inputGained_1) + (B2 * inputGained_2) - (a1 * output_1) - (a2 * output_2);

    return output;

And this is running at 44.1  kHz: What equation could I used to calculate the group delay at a given frequency like say 82 hz?

Thanks.


texts: Is there a noise cancellation technique for over-lapping spectra?

texts:
I understand Wiener filter is used to reduce noise, Gaussian noise as well, once we give the noise characteristic. One may derive the filter coefficient in adaptive way as well But can it remove noise if the noise is in the same band as the desired signal(EDIT: without loosin in-band energy of the signal)? Is there an such filter which is applicable for in band noise?  As per my understanding there is none and is not possible.
Edit: The signal is not assumed to be stationary.

Any suggestion, comment please?


texts: Difference between Sub Sampling and Down Scaling of Images

texts:
I know that there are many ways to upscale (interpolate) an image using bilinear, bicubic, sinc... Somehow, these same algorithms can also be used to downscale an image. But when it comes to subsampling, I've come across only two methods: Gaussian-blur then discard some rows and columns, or block-averaging.

Is there such thing as bilinear, bicubic,... subsampling?

Instead of Gaussian-blurring then discarding or block-averaging, can/should I do Gaussian-blurring then block-averaging? or bilinear, bicubic,...subsampling if there is such thing?

What I'm doing is trying to make image Gaussian and Laplacian pyramids; but I'm not sure how to decimate and interpolate the image. (Some Python code with Scipy/Numpy would really help.)

Thank you very much.


texts: Interpretation of a periodogram: the effect of window sizes

texts:
I am specifically working on a categorical time series of best n number of parallel threads (for n=1..16) for a multicore workstation. I constructed the Walsh-Fourier periodogram for size 2^10 and 2^11 data points. Both of the images are with the first value removed to allow a scaled image. The difference is dramatic as you can see.

My naive interpretation is that there is something stationary going on with the $2^{11}$ range that isn't happening in the $2^{10}$ range. In particular my data covers minute intervals and if I am looking for background programs that will effect the the best thread policy, if they are scheduled nightly (antivirus, back-up etc) then in terms of 1440 minutes=1 day they will fall within $2^{11}$ minutes because $2^{10}<1440<2^{11}$. Am I making the correct assumption here?

Because I am too new on this particular stackexchange I am not allowed to post the periodograms but they are hosted here.


texts: Transform one-pole cascaded system into a biquadratic section

texts:
I have two cascaded one-pole sections with the same coefficient b for both sections.
$$
y_A[n]=x_A[n]+b\cdot(y_A[n−1]−x_A[n])\\
y_B[n]=y_A[n]+b\cdot(y_B[n−1]−y_A[n])
$$

The final output of the system is $y_B[n]$.

I would like to re-implement the system as one biquad section of the form:

$$
y[n]=a_0x[n]+a_1x[n-1]+a_2x[n-2]-b_1y[n-1]-b_2y[n-2]
$$

The final output of the system is $y[n]$.

I would like to know what calculations I have to perform in order to make the output of the biquad exactly match the output of the one-pole cascaded system.
$$
y_B[n] = y[n]
$$
Obviously the 5 biquad coefficients $a_0\; a_1\; a_2\; b_1\; b_2$ will have to be determined as a function of the $b$ coefficient as found in the cascaded system above.

Can someone lead me through the math?

Thank you.


texts:
I'll start off by saying, yes this is a homework problem, and no I do not want the answer. I simply need to know where to start.

The input to a system with a frequency response of $H(\omega)$ is $x(n)=A\sin(\omega_0 n+\theta)$

What is the expression for the systems output, $y(n)$.

As I said, I just am unsure of where to start.


texts: DFT conjugate of $X^*[k]$, how to prove its formula in terms of $x^*[n]$?

texts:
Trying to prove that:

$$
X^{*}[k] = \sum^{N-1}_{n=0} x^{*}\left((N-n)\right)_N\ W_N^{nk}
$$

Where:
$$((x))_N \text{ = x  modulus N}$$

$$W_{N}^{nk} = e^{-j\ 2\pi / N}$$

So I start out with definition for DFT:

$$
X[k]=\sum^{N-1}_{n=0} x[n]\ W_{N}^{nk}
$$

Then I conjugate both sides of equation.

$$
X^{*}[k]=\left[\sum^{N-1}_{n=0} x[n]\ W_{N}^{nk}\right]^{*}
$$

$$
X^{*}[k]=\sum^{N-1}_{n=0} x^{*}[n]\ W_{N}^{-nk}
$$

then somehow the -n changes to N-n which i'm guessing is this:

$$((-n))_N = N - n$$

I don't really understand it... but, seems its equivalent if you do this inside of an exponent… to just randomly take a modulus of a negative signed number and replace it with its modulus equivalent?

$$
X^{*}[k]=\sum^{N-1}_{n=0} x^{*}[n]\ W_{N}^{(N-n)k}
$$

Then, I get stuck on the proof because somehow $x^{*}[n]$ changes to $x^*((N-n))$,

and $W_{N}^{(N-n)k}$ somehow changes to $W_{N}^{nk}$.

and thus, i'm unable to obtain final result of:

$$
X^{*}[k] = \sum^{N-1}_{n=0} x^{*}((N-n))_N\ W_N^{nk}
$$

still wondering what the rules are to apply to get there?


texts: Hilbert Huang Transform: Negative value in instantaneous frequency

texts:
I'm working on Hilbert Huang Tranform and have a question on computation of instantaneous frequency. I'm using HHT toolbox by Patrick Flandrin for EMD and computing IF manually. Problem is that I am getting negative values in IF.

Signal used: a $10\textrm{ Hz}$ sinusoid

 Ts=1/Fs;
 t=0:Ts:4;
 z= 5.*sin(2*10*pi*t);

Since this is a monocomponent frequency, IMF is the signal itself. I compute analytic signal using MATLAB function hilbert and then apply the following code to compute IF:

function [f] = inst_freq_manual(analytic_signal,Fs)

factor =  Fs/(2*pi);
N= length(analytic_signal);

% Differentiation formula
fin = factor*(angle(analytic_signal(2:N).*conj(analytic_signal(1:N-1))));

% Insert leading 0 in return-vector to maintain size
f = [0 fin];
end

I find that the IF values are negative. Mean of IF values =-9.9975 (correct numerical value but negative sign).

Plot of instantaneous parameters:

As per literature, the IMFs from HHT should have positive instantaneous frequencies.

Can anyone please tell me where I am going wrong?
Is using absolute value an option?


texts: Filtering an image with a helmholtz-type equation

texts:
I'm attempting to apply the following PDE as an image filter to smooth a discrete heightmap with a helmholtz-type equation as described in this paper. It seemed like an interesting alternative to a gaussian filter.

The equation is:
$$ddx(h') + ddy(h') + y(h'-h) = 0$$

I solved for h and got:
$$\dfrac{ddx(h')}{y} + \dfrac{ddy(h')}{y} + h' = h$$

I then discretized it with a central finite difference and turned it into a linear system of the form $Ax=b$ where $b$ is the source image and A is a matrix with elements around the diagonal corresponding to the coefficients of a central finite difference approximation of the second derivative. I also added an additional $1$ to the diagonal to account for the standalone $h'$ on the left hand side.

Unfortunately, the results don't look anything like a smoothed version of the original image. For very high y, the resulting image is mostly similar to h, but for lower y it quickly degrades into noise and eventually just a black image.

I suspect part of the issue is the way I'm dealing with the boundary conditions. I recognize that the derivative is undefined at the border of the image, so I've tried a number of different approaches to address this from excluding the borders from the system, to special casing the kernel of border pixels by adding the missing border weights to the diagonal of A.

Here is some code I've been playing around with to solve this problem.

I would greatly appreciate help understanding how to handle the boundary conditions and learning how to properly apply this PDE as an image filter!

Thanks,
Kris


texts: calculate SNR of two signals with power spectral density

texts:
You can see my system's diagram as shown below

At output,two different input signals are sited in each ($\frac{H}{3}$) position $(w_2=w_1+ \frac{H}{3})$ and if output power spectral density equation is included different sentences such as

Dc term reaction of input noise to itself.
reaction of inputsignal(w1 or w2) to itself.
reaction of input signals(w1 and w2) to themselves
reaction of input signal(w1 or w2) and input noise
reaction of input signals(w1 and w2) and input noise to themselves.

How can obtain input and output SNR?


texts: Most common modern day Image Segmentation techniques

texts:
I was doing some reading up on image segmentation techniques and I was wondering about modern day, state-of-the-art segmentation algorithms.

What current segmentation techniques that are a 'must-read', i.e. currently most commonly used in the community? Which techniques did you came in contact with and found most effective and useful (and for which application)?


texts: I am working to understand and use the Chirp Z-Transform.  I want to use the algorithm for simple signal processing on data sets that are not a power of two.  I need to be able to inverse transform as I want to transform a set of data to the frequency domain and operate on the (complex) frequency coefficients, and then transform back to time.  So for example, one operation I perform is to integrate the signal using the Omega method.  Once I operate on the coefficients I need to transform back to the time domain.  I believe I have the Chirp Z-Transform algorithm figured out.  I can for example, pass in a 100 point time series made up of a few tones, and transform the series to the frequency domain using the Chirp Z-Transform.  I get spectral content right where it should be.  The question is, what do I do with those 100 complex valued coefficients to transform them back to a time series.  At first I thought I would use a power of 2 FFT with the additional coefficients padded to zero.  That didn't work and I now understand why that isn't the same thing as transforming only the 100 coefficients back to time.  Could someone help me to understand how to inverse transform a non-power of 2 set of Fourier coefficients.


texts: Fourier Spectra : Significance of the Negative Amplitude

texts: For example, for an aperiodic gate pulse, the Fourier Transforms for the continuous time case is a sinc function, while the discrete time case gives a sine over sine periodic kind of a function. In either of these cases, we get a spectra that goes below the frequency axis, i.e., acquires a "negative amplitude" (namesake) on many occasions in its excursions. I realize that the spectra gives a measure of "how much" a frequency component is present. But how this negative amplitude satisfies this criteria is what I am missing. Is it kinda like.. "How much less that frequency component is present"..? I know I might not be making any sense. Please help. Thanks!!


texts: Pyramid vs Scale-space in Focus Stacking

texts:
I'm working on a focus stacking program, I wrote one a few years ago which used a very simple method, basically taking a DoG approach, but only operating on the single full resolution image. The resulting images had a soft focus look, kind of dreamy.

This time I researched further and found papers describing image pyramids, essentially doing the same thing as before except on multiple spatial frequencies, then recombining them with the maximum components which I imagine will solve the soft focus look.

But learning about pyramid based approaches led me to scale-space representation and multi-resolution analysis, and I can't tell if I should pursue one of those methods for possibly better results? (And MRA led me to wavelet transforms which I also don't see how I would implement.)

Furthermore, I don't really understand how either would even be implemented for focus stacking, but my impression of scale-space is that in practice I would just end up building a pyramid anyway, but that doesn't seem right. So am I missing something about that?

And does it seem like there is anything to be gained by using a newer or more sophisticated approach than image pyramids for focus stacking?

(I used Python last time but am using Apple's CoreImage framework this time, so I might be a bit more limited now, maybe MRA isn't even an option.)

Thanks in advance!


texts: How the Grayscale images are converted back to RGB?

texts: In case of RGB to grayscale conversion we take the average of R,G,B channels and convert them to Grayscale which has songle channel. How to convert back the Grayscale to RGB ??


texts: How are the constellation point coordinates determined in digital modulation

texts:
I need to understand where the concept of the distance comes from. What I mean is, how do we know what coordinates are to assign to points in the constellation. For example here for QPSK the points are given coordinates $$\pm\sqrt(E_b),\pm\sqrt(E_b)$$
For BPSK on the other hand they are $$\pm\sqrt(E_b),0$$
This question was answered here as:
"A symbol with coordinates $(x,y)$ has energy $x^2+y^2$", and then it makes sense as for example for BPSK the energy per symbol is just $E_b$ and $\sqrt(E_b)^2+0^2=E_b$ but where does this come from? How to derive this? Thanks.


texts: Can you explain likelihood function in the equation?

texts: In the paper here, just before equation 12, authors introduce likelihood function $\pi(r|\theta)$. What does $\theta$ represent here ?. What is the intuitive explanation behind equation 12?


texts: Possibility of using Power allocation with Maximum likelihood estimation

texts:
I have a basic question about maximum likelihood (ML) estimator and its implementation with power allocation.
Suppose we are using $MIMO$ system with $4$x$4$ transmit and receive antennas, so the received signal which includes the signal + noise + ISI can be written :
$y = HX + n$
$H,X, n $ are the channel, transmitted signal and noise, respectively.
suppose that only first antennas is active and other antennas are not active, so the transmitted signal will be as below:

My question is about the possibility to use power allocation to detect the active antenna, which means the estimated signal will be as:

I am asking if we can detect the active antenna in that way, is that possible ? Is there another method we can use instead of that ? Is there any document explain that ?
Thank you.


texts:
Let's assume I send a root raised cosine pulse of 1ns at t=0 through a wireless channel.
The power measured @ Receiver from t=0 is shown in image.
Now from this image, I know that there are multipath components(due to spreading of signal power across time).

I want to know that up to which time do I have to measure the power for a single symbol transmitted and how can I know that after a certain amount of time the second symbol will start. (how can I measure that certain amount of time?)
In this received power I know that(by spreading) that there are multipath components. I want to measure the Gain of those multipath components. So how can I sample it and how do I know about the number of multipath components?


texts: Cross Power Spectral Density of Unevenly Sampled Data

texts:
Here's my problem. The input signals $x$ and $y$ will be having the time value aligned with each other. However, the data are not evenly sampled. I would like to calculate CPSD of both signals.

The solution comes to my mind is as follow

$R_{xy}$ = cross correlation of $x$ and $y$ (I'm not sure how to do it with unevenly sampled data)
Use Lomb Periodogram to estimate the CPSD of Rxy

Is there anything I miss? And what is the right way to do it?

Update

Sam Maloney suggested to use interpolation to fill the gap to produce evenly sampled data.
This solution is good.

However, the data I obtained may have some missing data / long gap. Interpolating the data may "predict" the data wrongly and produce undesired end result.

As quoted in Numerical Recipes (Spectral Analysis of Unevenly Sampled
Data)

  However, the experience of practitioners of such interpolation techniques is not reassuring. Generally speaking, such techniques perform poorly. Long gaps in the data, for example,
  often produce a spurious bulge of power at low frequencies (wavelengths comparable to gaps)

This is the reason why I choose Lomb Periodogram over FFT in Step 2 without interpolating the data. However, I'm unable to figure out how to solve the Step 1. Is there any way to calculate the cross correlation of 2 signals without interpolating them?


texts:
MCQ
The region of convergence (ROC) for a right-sided (infinite) sequence is
1 point
Entire z-plane excluding zero and infinity
Interior of a circle
Entire z-plane including zero and infinity
Exterior of a circle


texts: What does the frequency axis of a Power Spectral Density mean?

texts:
I have never really understood what the frequency axis meant when we plot the Power Spectral Density(PSD).

Does it correspond to frequency as we get after we take the Fourier Transform of a time domain signal. This doesn't make sense to me because the frequency here(when we do FT of a time domain signal) corresponds to the time variable "t".

To find the PSD we take the Fourier transform of the Autocorrelation function which is a function of either time difference $$\tau = t_2 - t_1$$ i.e. $$R_{xx}(\tau) $$ (for the Wide Sense Stationary case) or $$R_{xx}(t_1,t_2)$$ for other cases.

Here the autocorrelation is a function of either time difference or two time variables. If this is the case then how can we say that the PSD gives the power distributed in the frequency domain?

How do I understand this. If someone could give both a mathematical explanation and an intuitive explanation I would be really grateful.

Thank you.


texts: Applying Photoshop's "Shadow / Highlight" Correction Using Standard Image Processing Algorithms

texts:
Adobe's Photoshop has an excelent filter called "Shadow/Highlight" used to correct, in one step, dark areas AND overexposed areas of a photo. By simple analysis, I suppose it uses some radius-dependent adaptive contrast algorithm (halo effect is typical), possibly combined with interpolation (even extra-dark or extra/white, clipped areas get very improved with this filter.

The parameters used can be seen in the dialog box below

What I would be looking for is a way to start exclusively with the parameter values given in the dialog box, apply a sequence of parameterized operations "blindly", and get back the result non-interactively (well, perhaps previewing the image while adjusting parameter sliders in some GUI environment).

Any suggestion about how to perform this with standard, well-known image-processing algorithms?


texts: How does sampling rate of $x[n]$ relate to sampling rate in frequency domain after DFT transformation?

texts:
I've got an analog signal $x(t)$ sampled at frequency $F_s$ to obtain samples:

$$ x[n] = x(t) \bigg|_{t=n/F_s} $$

I transform this signal with DFT defined as:

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-i2\pi kn/N}$$

Is there such a thing as a sampling frequency / rate in the frequency domain? How does it relate to $F_s$?


texts: What is the difference between i.i.d noise and white noise?

texts:
I want to know the difference between independent and identically distributed (i.i.d) noise and white noise.
In my short knowledge, i.i.d is that there is no relationship about time dependancy. White noise means that there are relationship about time dependancy.
Actually, I'm not sure whether this is correct or not. Also I want to know  what is an iid white noise.
Can you tell me where we find the iid noise in the nature?
update
enter image description here
enter image description here


texts: Demodulating signal by passing incoming carrier wave through low pass filter as opposed to multiplying by carrier wave?

texts:
Demodulation through many resources found online says that the incoming signal is multiplied by the carrier wave and then passed through a low pass filter that removes the carrier wave from the modulated signal. I understand this.
However, we know that the incoming signal is already multiplied by the carrier wave therefore, shouldn't passing the incoming signal through a low pass filter remove the carrier wave anyway?
For example.
The former states :
Modulated Signal * Carrier wave  = incoming signal
Incoming signal * (Wave of same frequency and phase)  ---> low pass filter = modulated signal
The latter states(my question):
Modulated Signal * Carrier wave  = incoming signal
Incoming signal ---> low pass filter (removes the carrier wave)= modulated signal
I would like to know why the former is stated rather than the latter and the benefits of it.
I suppose the incoming signal can include noise that could somewhat alter the signal therefore multiplying by a wave of same frequency can help, but I'm not too sure about this.


texts: Is this how the lowpass butterworth analog filter can be implemented?


texts: What are some face detection (not recognition) algorithms suitable for limited (embedded) processors?

texts:
I am looking at a hardware based platform for surveillance class of applications. I want to identify faces of people as they appear in the scene.

I don't intend to really apply face recognition (then or later). Also orientation of faces can be assumed straight for practical purposes. Face detected once needs no further processing.

Is there is a fast (efficient) algorithm to detect faces from the live video source?


texts: How to calculate the phase shift AND time delay of non-periodic signals

texts:
Say I have two non-periodic signals, $f_1(t)$ and $f_2(t)$, with Fourier transforms $F_1(\omega)$ and $F_2(\omega)$. Basically, I need to line up $f_1(t)$ and $f_2(t)$ as close as possible, and I am allowed to shift the signals in time and multiply their Fourier transforms by a constant phase.

What is the most efficient method to find the values of the phase shift $\Delta \phi$ and the time delay $\Delta t$ such that, when the Fourier transform of $f_2$ is modified into

$$ F_2(\omega) e^{i (\Delta \phi + \omega \Delta t)}, $$

the modified $f_2$ best "lines up" with $f_1$ according to their cross-correlation (or some other measurement)?

I am interested in simple and computationally efficient methods, even if they are not necessarily perfect.


texts: Modelling propagation of sound wave by particle simulation

texts:
How might one model the propagation of a sound wave through a 3D environment with obstacles in it?

My initial thought would be to use freely moving particles that repel one another; with sufficiently many particles I see no reason why the simulation wouldn't work.

However, if any particle could collide with any other then this introduces a huge complexity hit.

So I'm wondering whether the particles could be effectively anchored on a sphere packed lattice (dodecahedron?), so each particle has 18 neighbours. (you could think of placing one tennis ball on the table, then six fitting around it, now 6 can be fitted on top, and by symmetry 6 could be fitted underneath. so the original ball has 18 neighbours.

And then shockwaves could propagate through this mesh, and at each iteration each particle just needs to calculate the resultant force of its 18 neighbours acting on it.

I can't see whether this model would work or not.  I suspect it is going to model sound propagating through a diamond rather than sound propagating through air, but isn't it going to be essentially the same?

Also, this isn't quite on topic for DSP, my apologies. Can anyone recommend a more suitable place to ask this question?


texts: When calculating the Nyquist frequency should carrier frequency be included

texts:
This isn't my field, so please forgive any misnomers on my part.

I'm looking into buying a amateur radio receiver, and I'm thinking about what sampling rate I need. I'm most interested in satellite signals, which transmit on ~400MHz (for example). So that would imply I have a nyquist frequency of ~800MHz, or 800 million samples per second. This seems much higher than any hardware I can find.

My current thinking is that the frequency to be used in the nyquist calculation isn't the carrier frequency. Which frequency should I use?


texts: What kind of edge detection use paint.net?

texts:
What kind of edge detection is used Paint.NET program? It allows me to choose the angle of the edge not like the one used by opencv (Canny).

Here an example:


texts:
I was calculating which would be faster in computation, Inverse Fast Fourier Transform or Histogram Of Gradient but I am unable to calculate its running time.

What would be the running time of Histogram Of Gradient if the image has 'p' pixels ?


texts:
I come across different types of bandwidth. In digital modulation should it be null to null always?
when do we use 3 db bandwidth?


texts: This is the first time I'm diving a little deep into images and computer vision. I have an aerial image taken from a drone (I know the pitch, roll, yaw of the drone and the angle of the camera relative to north). How do I transform the image so that it looks as if it is taken from the top? I want to map latitude, longitude on the aerial image to pixels (I know the lat, lon of all four corners of the image) and my understanding is that I need to transform the image before I can interpolate and map a new lat,lng pair to pixels.


texts:
Analogue radio stations often have a range of frequencies that they can use. Different transmission towers will transmit on different carrier frequencies in this range, because the signal takes time to propagate through air and if they used the same carrier then the received signal would be... Well, it wouldn't be well-formed FM, and it would almost certainly be echoey and distorted.

But what if they did all use the same carrier? You'd need multiple aerials placed some distance apart (would 3 be enough?) and would have to do something clever with them... But would it be possible to extract a signal from that?

If it's not too different a question, would it be possible to extract the individual signals if the transmission towers were transmitting something different to each other on the same carrier frequency?


texts:
I am new to the study of time series. Recently I have asked a question about the covariance of real and imaginary part of a real(in time domain) stochastic time series  and I have received an answer for it. The problem is that for continuous time series the variance of each point in frequency domain is infinite. I've been told there that its the reason that they use power spectrum. Now the confusion comes from the fact that I don't know whether power spectrum is time averaged or not? More precisely as wiki says:
$$S_{xx}(\omega)=\lim\limits_{T\to \infty}\mathbf{E} \left[ | \hat{x}_T(\omega) |^2 \right] = \lim\limits_{T\to \infty}\mathbf{E} \left[ \frac{1}{T} \int\limits_0^T x^*(t) e^{i\omega t}\, dt \int\limits_0^T x(t') e^{-i\omega t'}\, dt' \right] = \lim\limits_{T\to \infty}\frac{1}{T} \int\limits_0^T \int\limits_0^T \mathbf{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt$$
OR
$$S_{xx}(\omega)=\mathbf{E} \left[ | \hat{x}_T(\omega) |^2 \right] = \mathbf{E} \left[  \int\limits_0^\infty x^*(t) e^{i\omega t}\, dt \int\limits_0^\infty x(t') e^{-i\omega t'}\, dt' \right] =  \int\limits_0^\infty \int\limits_0^\infty \mathbf{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt$$

For example in Wiener-Khintchine theorem as far as I can see there is no time averaging:
$$r_{xx} (\tau) = \int_{-\infty}^\infty S_{xx}(f) e^{2\pi i\tau f} df$$

And is there any difference when the signal is discrete?


texts: What is a correct approach to baud-rate digital timing recovery for self-equalized line code

texts:
In particular, I have a signal coming from a twisted pair of arbitrary length with a differential encoded
biphase line code. What should be a correct approach to timing recovery in this case?
Currently I am trying to implement a baud-rate timing recovery algorithm using a Mueller&Muller TED (timing error detector), loop filter (average over 120 symbols) and counter-divider with adjustable initial value as a digital VCO. The source code for current implementation of the algorithm in python is here. For TED I use the following timing function:
$$
\frac{(a_k^2 - 5)}{16}(a_k x_{k-1} - a_{k-1} x_k )
$$
For simulation purposes I have sampled versions of the signal in half-duplex mode (i.e. without echo). Signal samples look like this for two different channel lengths:

For these signals the timing recovery algorithm converges to the phase near signal's zero crossings. This does not allow to correctly determine the original bit sequence (e.g. "+1 -1 +1 +1 -1 +1 ...") and I would like it to converge to the maximum eye opening, i.e. somewhere near the peaks.
I've managed to get this by adding a $(1 - aD)$ filter, where $a$ is $0,75$ for the first signal sample and $1,325$ for the second one. But because a channel has an arbitrary length, the shape of a signal is unknown during design phase and this filter with fixed coefficients seems like a bad solution. Moreover, adding this filter was purely an empirical solution and is not backed up by any math.
I am very new to timing recovery and would appreciate any advice.


texts: How to plot poles, zeros and impulse response of a transfer function

texts:
I want to plot poles and zeros as well as impulse response of a transfer function that is given below:
$$H(z) = -0.0625z^4 + 0.25z^3 + 0.625z^2 + 0.25z -\frac{0.0625}{z^4}$$
For poles and zeros, I used MATLAB function residuez. However, I got error which is similar in impulse response plot.
code:
 b = [-0.0625 0.25 0.625 0.25 -0.0625];
 a = [0 0 0 0 1];
 [R,p,C] = residuez(b,a); % for poles
 y = impz(b,a); %for impulse response
 figure(1);
 stem(y);

Could someone help.
Thanks


texts: Calculate Frequency/Period from Output of a Maximum Entropy Spectrum MATLAB function

texts:
When I plot said function output, "p," it looks like a typical spectrum with an easily identifiable peak at a given ix on the x axis. The help rubric from the function is

SYNTAX: p = mesa(x,m,nfreq);
For a vector x, this function calculates a maximum-entropy spectrum
of order m. The spectral estimate is returned in the vector p, which
has nfreq points linearly spaced in the Nyquist frequency interval 0-.5.
The psd is normalized such that the mean square value of x equals the
integral of p from -.5 to .5, so sum(x.^2)/N ~= sum(p)/nfreq.
Mesa is based on the Burg algorithm, as described in Numerical Recipes
and implemented in their memcof and evlmem subroutines.

Written by Eric Breitenberger        Version 5/24/95

My problem is that the x axis is unitless, i.e. if nfreq = 200 then the x axis is 1 to 200 long. If the peak is at x axis ix = 22, for example, how do I calculate the measured frequency/period of the input signal "x"?

To give a concrete example:

the length of the input vector "x" is 20 because I believe (or know due to a controlled input experiment) this length vector contains one full cycle (noisy sinusoidal) of interest

order m is set to 10, for example, although this is incidental to my question as I can see from experiment that changing the value of m only affects the sharpness of the peak

nfreq is set to 200

and with these values the peak occurs at x axis ix == 22.

From this, what is the mathematics to get the measured frequency/period from the output of the function?

Plot of the output

for a known input, e.g.

input_x = awgn( sinewave( 20 , 20 ) , 10 ) ;


texts:
Let's say I have 2 sequences a and b in the time domain. Both are length N. A and B are the DFT of a and b. If I do a circular convolution of A and B in freq domain (A o B), then the IDFT of the result will be a x b in the time domain. Correct?
Now what happens if I do linear convolution of A and B in the freq domain, i.e. A*B? The results will have a length of 2N-1. What does the IDFT of the result mean in the time domain?
Notation:
-x is multiplication
-* is linear convolution
-o is circular convolution


texts: Some problem with the expression of QAM basis and signal

texts:
The expression of QAM basis is
\begin{align}
\psi_1&= \sqrt{\frac{2}{E_g}}g(t)\cos(2\pi ft)\\
\psi_2&= -\sqrt{\frac{2}{E_g}}g(t)\sin(2\pi ft)
\end{align}

I want to ask why is $\psi_2= -\sqrt{\frac{2}{E_g}}g(t)\sin(2\pi ft)$,
instead of $\sqrt{\frac{2}{E_g}}g(t)\sin(2\pi ft)$? Why is there a "-" (i.e. minus) before $\sqrt{\frac{2}{E_g}}$?
And QAM signal is $\ s_m(t)=A_{mi}\sqrt{\frac{E_g}{2}}\psi_1 +A_{mq}\sqrt{\frac{E_g}{2}}\psi_2$, does the $\sqrt{\frac{E_g}{2}}$ here have a relation with its basis $\sqrt{\frac{2}{E_g}}$ ?


texts:
What fading model should I use for simulating rain loss in my network simulations? I am using NS-3 for this and I'm quite lost on how to go about this.

My simulation is just simple packet streaming and I want to test the effects of various interference on throughput.


texts: I'm seeing everywhere that the DFT of a sinusoid (i.e. a sine that has been windowed with a rectangular window) are just two deltas. What I don't understand is why the DFT is not the convolution of the transform of the infinite sinusoid (two deltas) and the transform of the rectangular window (a sinc). Wouldn't this be the right thing to do?