texts: Accelerometer - coordinate system transformation

texts:
I'm getting some accelerometer readings from an Android phone, but it comes in on the phone's coordinate system. I want to apply a transformation to put the acceleration in the world coordinate system (e.g. so Y will always be up and down, regardless of how the phone is rotated)

Android allows you to get the accelerometer values. It can also calculate a rotation matrix (using a gyroscope and a magnetic sensor), so all I need to do is multiple the acceleration vector through by this matrix and I have transformed values. I followed the steps here: https://stackoverflow.com/questions/14963190/calculate-acceleration-in-reference-to-true-north

However, there are 3 sensors at play here: acceleration, gyro, magnetic. The way android operates is event based, so whenever it detects a change in any of these sensors, it will give you back a value.

Given these are 3 different sensors, the events for each one is very likely to come in at different times, and you may get more readings from one sensor than another within some time period

My question is, does varying amounts of spacing between readings for the 3 different sensors matter for doing coordinate system transformations? And if so, what is the best solution for doing this type of transformation with event based sensors?


texts: Compare phase and magnitude spectrum results of 2 Images

texts:
I have read that from Fourier transform we obtain magnitude and phase spectrum. The magnitude spectrum tells you how strong are the harmonics in a image and the phase spectrum tells where this harmonic lies in space.

I have used Matlab to compute the magnitude and phase spectrum of the 2 grayscale images of babies using fftn function, but the all spectrum are difficult to understand.

Can anyone please explain the features given in 1st paragraph by comparing the spectrums of both babies ?

1.phase spectrum

and

2.magnitude spectrum

(If you have results other than in Matlab, explaining with your method is also okay)

image1=imread('D:\baby1.jpg');%Read colour image
image1=rgb2gray(image1);%covert colour image into Gray scale

image2=imread('D:\baby2.jpg');%Read colour image
image2=rgb2gray(image2);%covert colour image into Gray scale

figure,
subplot(1,3,1);
imshow(image1);
title('Gray scale Image of Baby1');

%Apply Fourier transform on Gray scale image
fft1=fftn(image1);
fft1=log(1+fftshift(fft1));%Use log function for scaling

%Find the magnitude spectrum
magnitude1=abs(fft1);
subplot(1,3,2);
imshow(magnitude1,[]);
title('Amplitude Spectrum 1');

%Find the phase spectrum
phase1=angle(fft1);
subplot(1,3,3);
imshow(phase1,[]);
title('Phase Spectrum 1');

figure,
subplot(1,3,1);
imshow(image2);
title('Gray scale Image of Baby2');

%Apply Fourier transform on Gray scale image
fft2=fftn(image2);
fft2=log(1+fftshift(fft2));%Use log function for scaling

%Find the magnitude spectrum
magnitude2=abs(fft2);
subplot(1,3,2);
imshow(magnitude2,[]);
title('Amplitude Spectrum 2');

%Find the phase spectrum
phase2=angle(fft2);
subplot(1,3,3);
imshow(phase2,[]);
title('Phase Spectrum 2');


texts: Steady-State Output from Transfer Function

texts:
The progress I have made is as follows:

$\sin(t)$ is our signal therefore $\omega = 1 = 2\pi f$ and $f$ = $\frac{1}{2 \pi}$

Also, $f_s$ = 10Hz therefore T = $\frac{1}{f_s}$ = 0.1s

$H(z) = \frac{z}{z-0.7}$ so let $z = e^{j\omega T}$ where $\omega T = 1 \times0.1 = 0.1$

$H(\omega) = \frac{e^{j0.1}}{e^{j0.1} - 0.7}$

Taking the denominator:

$\cos(0.1) + j\sin(0.1) = 0.995 + 0.099j$

therefore,

$\cos(0.1) -0.7 + j\sin(0.1) = 0.311∠0.324$

$H(\omega) = \frac{1∠0.1}{0.311∠0.324}$

$H(\omega) = 3.21∠0.224$

From here I am out of ideas on how to continue.  Any advice appreciated.


texts: Determining Stability of a continuous time system using Laplace Transform

texts:
I'm following Oppenheim's book. In exapmles, the laplace transforms of of the following signals

$e^{-t}u(t)$ and $e^{-t-1}u(t+1)$

is given as $\frac{s}{(s+1)}$ and $\frac{e^{-s}}{(s+1)}$ both having the same ROC $Re(s)>-1$.

however, the first signal is causal and the second one is non-causal. The system function is not rational for the second one.But why is that the problem when the ROC is essentially same?


texts: coded filter vs fdatool : inconsistency in the results for same parameters

texts:
I tried creating a filter in 2 different ways in Matlab and I obtain different results: one filter is stable, whereas the other one is unstable.

Coded filter

delta_t = 1.53846e-04;
Fs = 1/delta_t;

cut_F = 8;
Wn = cut_F/(Fs/2);
ftype = 'high';
[b,a] = butter(8,Wn,ftype);
fvtool(b,a)

The result is an unstable filter

fdatool

This is what I obtain with the toolbox using the same parameters:

I am using already implemented functions in both cases...Why I get so different results?


texts: Is it possible to design a band-pass filter with fewer taps than a series HPF + LPF?

texts:
I am digitally filtering a signal that is streaming live. The goal is to reduce the filter delay as much as possible while maintaining the quality of some measurements based on the filtered signal. For legacy reasons, the system was built with a low-pass and high-pass FIR filter in series - each filter has the same number of taps (kernel size), $n$.

I realized though, we are no longer constrained to these and could implement a single FIR band-pass filter. If I can reduce number of filter taps involved from $2n$ that will reduce the delay. My question is: is it possible to design this band-pass filter with fewer taps (a smaller kernel) than the combined taps (kernels) of the LPF/HPFs while achieving the same performance?

My initial thought is FIR filtering is a convolution operation, and thus linear + associative. That means you can construct the band pass kernel by convolving the LPF + HPF kernels. I'm unsure how to interpret this in terms of taps/kernel size because the convolution in some sense assumes an infinite kernel.

Thanks!


texts: Are there paper about black box testing of optic flow algorithms?

texts:
Do papers exist about comparisons of the resulting optic flow and its precision compared to the used optic flow method?

I would like to find paper that dealt with this topic in some way and had discussed different ways to conquer this problem.

Thank you in advance.


texts: Finding correct values based on information from two arrays

texts:
Consider the following scenario:

Say, one machine is sending out a beep signal every 10 seconds in a very noisy environment. I have two sensors which detects these beeps independently. Device A is near to machine compared to B, so B always have a two samples delay compared to A. Since it is noisy environment, devices sometimes detects noise also as beeps and sometimes even misses the data.

Below is the data

A = [5,10,20,25,30,40,45,50]
B = [12,18,22,32,38,42,52]

See, A and B have beeps at 10 & 12 which says beep @5 in A is a noise. Similarly @25 in A because 30 in A and 32 in B says they are the good data.

Similarly in B, 18 is a noise because there is a synchronized peak at @20 in A and @22 in B.

Below is a plot of that (red mark shows noisy beeps):

Question:

Is there any mathematical/signal processing method which I can use to find the location of correct beeps based on the synchronized information available from two datas and get rid of the noisy ones?


texts:
I am trying to produce some simple audio waves using WASAPI (Windows Audio Session API). If I skip all the abstractions, my program basically receives a buffer from the API, then writes results of a sine function (or some other simple function) of a selected frequency to the buffer sample by sample, while counting the samples to keep track of time. When it goes to the next buffer it continues with an offset of the samples passed. So in theory, it should be a simple quantization of the sine wave.

The problem is, that some frequencies work and some don't. In particular, multiples of 50 work perfectly fine, while all the other frequencies sound glitchy.

Any ideas on what problem this could be related to? Is it some quantization error maybe? Should I look into some concepts like signal-to-noise ratio? Or are these problems a matter of the DAC and not my concern.

Edit:

I agree with the comments bellow, that my code being wrong is more probable, but after trying to inspect every detail of the code and also noticing this strange behavior I just thought I might ask some more experienced people, and hope that this behavior might be suspicious of some more general problem. If it sounds like a problem specific to my implementation, then I'll try to post a simplified example of my code. (Only if this community is a right place to post such questions of course.)


texts: What is a correct way to downsample an impulse train

texts: So I have a signal which is caused by an impulse train. I can determine the impulse train with 500hz but sample my output signal with only 10hz. Then I want to calculate an impulse response function. So I have a binary "signal" for the impulses with 50 times more samplepoints then my actual signal. If I resample my binary signal using fourier which I would do with an ordinary signal, I get Gibbs ringing and I doubt my impulse response function would still make sense. I can choose the nearest neighbor to be 1 in my "downsampled" impulse train but then the timing is less accurate, I could divide the impulse over the two neirest neighbors bu then no single point will have the value of 1. What would be a correct way to go about this?


texts: Plot Spectrogram and FFT plots from samples collected by Atheros AR9271

texts:
I'm entering the world of signal processing, and as one of my first projects I've decided to plot a spectrogram and FFT Plot.
The samples comes from an Atheros AR9271. More details can be found here: The AR9280 and later NICs have the ability to report FFT data from the baseband
Currently, I have a small script, similar to this.
What I wanted to do, is to use the raw abs(|I|+|Q|) which is given for each of the 56 bins and build my data from there.
According to 1 I can do the following:

bin_sum = 10*log(sum[i=1..56](b(i)^2)
power(i) = noise_floor + RSSI + 10*log(b(i)^2) - bin_sum

So I have power[dBm] for each of the 56 bins.
From here I simply append the power series for the 56 bins to the overall samples series.
samples = np.append(samples, power)

This might be mistake #1, appending the samples as such.
I have tried to calculate fft(using scipy):
center_frequency = 2432 * (1e6)  # Channel 5 center frequency
half_bandwidth = 20 * (1e6) / 2
fs = 2 * 2.4 * (1e9) # **Maybe mistake #2**
N = samples.size
T = 1.0 / fs
yf = fft(samples)
xf = np.linspace(center_frequency-half_bandwidth, 1.0 / (2.0 * T), N // 2)
plot(xf, 2.0 / N * np.abs(yf[:N // 2]))

And regarding the spectrogram(using matplotlib)
specgram(samples, Fs=1.0 / fs, cmap='plasma')

Problems:

I'm almost sure I'm making at least few mistakes and I don't know where.
The x-axis for the spectrogram does not match the frequencies I'm sampling(should it?). In the case of the FFT Plot I made the adjustments to xf.

Question:
What mistakes and bad assumptions am I doing?
Regarding what I've wrote, which I'm not sure is correct, here are my understandings

fs - Why I've used: 2 * 2.4 * (1e9)? Due to 802.11n on AR9271 and I'm working on the 20 mode, which is 2.4GHz with bandwidth of 20Mhz, I was thinking that nyquist theorem is applicable here, is it?.
xf adjustment I've made so that it will plot around the center frequency. I'm not sure if its the right way to do so.
And most importantly, is the samples array I construct as I did, which is simply appending the power array so that I have: samples = [bin0Power_fftFrame0, bin1Power_fftFrame0, ..., bin55Power_fftFrame0, ...... , bin0Power_fftFrameN, bin1Power_fftFrame0, ..., bin55Power_fftFrameN]

Any help over one of my misunderstandings/mistakes would be really helpful.


texts: Soft question: Why do we need to reconstruct a signal?

texts:
I am currently reviewing the Nyquist criterion which says I need to sample the signal at 2 times the maximum frequency in order to avoid aliasing in the reconstructed signal.

My question is that in all the applications I've done, I rarely ever remember having to reconstruct the signal after sampling it.

For example, in controlling a digital controller, all I have to do is to sample the signal and send the sampled signal to the controller. End of story, no reconstruction required.

In other applications, I sample the signal to analyze its frequency content. After looking at the frequency distribution, the signal is completely discarded. No reconstruction required.

Can someone make it clear why do we need to respect the Nyquist criterion if no reconstruction is needed (so no worry about aliasing effect)?


texts:
I have a series a of values (0 and 1) coming from a Brownian process with drift for which I am studying the autocorrelation.

I used two methods:

1) numpy autocorrelation:

corr = np.correlate(a,a,mode='full')/a.size
corr = corr[corr.size//2:]

2) Fourier transforms (Wiener-Khinchin):

A = np.fft.fft(a)
S = np.conj(A)*A/a.size
c_fourier = np.fft.ifft(S)

However, I get different results:

I am not very experienced with stochastics or signal processing, so I have some difficulties understanding where the difference comes from.


texts: How to find phase values on precise pixel locations of a given Image?

texts:
I have computed magnitude and phase spectrum of cameraman Gray scale image  using fft function. Here, we get magnitude and phase spectrum of the whole image. .Please see the reconstruction of image from phase only part from the accepted answer of Image Reconstruction:Phase vs Magnitude  which says  that in an image neighboring pixels of the edge,lines and corner are in phase . So  I want to find phase values on precise pixels

Given an image of dimensions 256*256, how can I find phase values of pixels at locations(100,121),(100,122) and (100,123) ?

Whether it is possible using Fourier Transform or not? If not ,what is the other way to do it? Can anybody explain it with (or without )  code?

clc;
clear all;
close all;

i=imread('C:\Users\RK\Desktop\cameraman.gif');
%i=rgb2gray(i);
i=uint8(i);

figure,
subplot(1,3,1);imshow(i);
title('Cameraman Gray scale Image');

f1=fft2(i);
f2=log(1+fftshift(f1));

m=abs(f2);
subplot(1,3,2);
imshow(m,[]);
title('Magnitude Spectrum');

phase=angle(f2);
subplot(1,3,3);a
imshow(phase,[]);
title('Phase Spectrum');


texts: Bit accuracy improvement versus sample averaging

texts:
Supposing we have multiple samples of a scaled signal $x$ given by $y_i=a_ix+n_i$ where $n_i$ is AWGN how does the bit accuracy of estimation of $x$ improve with averaging of samples $y_i$ (number of bits per sample)? Assume we know $a_i$.

Say I want $x$ with an error affecting only $m$ bits does $\log m$ samples of $y_i$ suffice?


texts:
Can anyone please suggest books or resources related to variational mode decomposition (with matlab application)? It was suggested to me to understand the intuition behind ADMM first. However, I am struggling with that too.
Also, I am kinda started to understand EMD, but if you have any suggestions, I will be happy to get them.
Thanks!


texts: What is the meaning of inverse signal bandwidth?

texts:
I have been studying the two ray model in communication systems.

The text I'm studying has two signals, $u(t)$, which is Line of Sight and $u(t-t')$, which is the reflected signal.

It quotes,

  If the transmitted signal is narrowband relative to the delay spread ($B \ll \frac{1}{t'}$) then $u(t)$ is approximately equal to $u(t-t')$.

What is the meaning of this statement?
How is the delay related to the bandwidth?

Sorry if the question wasn't clear the first time.
For the source, one can refer to https://en.wikipedia.org/wiki/Two-ray_ground-reflection_model or Andrea Goldsmith's Wireless Communication (mine is page 34, chapter 2) which quotes the above.

Inverse signal bandwidth is ($\frac{1}{B}$) and inverse delay spread is ($\frac{1}{t'}$). So either ($B \ll \frac{1}{t'}$) or equivalently ($t' \ll \frac{1}{B}$).

So let me rephrase. How can I show that $u(t)$ is approximately equal to $u(t-t')$ given either ($B \ll \frac{1}{t'}$) or equivalently ($t' \ll \frac{1}{B}$).  In other words, what is the proof of the sentence in quotes (or highlights)?


texts:
I am stuck in a small question which I am trying to solve.

The signal $f$ defined by $f(t) := \text{sinc}(t)$ is to be sampled. How should I draw the illustrations of the $T$-sampled signal $f_{T} = (\mathbb{Z} \ni n\longmapsto f(nT)) $ for $T \in \left\{1,\frac{1}{2},2 \right\} $?

Thanks.


texts: Recovering a Differential Equation From the Transfer Function of a Cascaded System

texts:
With respect to the below discussion, consider that we are talking about LTIC systems characterized by constant coefficient ODEs.

Consider a cascaded system whose transfer function H(s) is given by

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

"Recovering" the defining differential equation from this transfer function yields

$$(D^2 + 4D + 3)y = P(D)x$$

where $P(D)$ is unknown but irrelevant for the purposes of this discussion.
What are the characteristic modes of the system?

Is it acceptable to simply say that the natural response $y_n$ is given as below?

$$y_n = Ae^{-t} + Be^{-3t}$$

It would seem to me that the answer is no. Why can we recover an ODE from a transfer function in general. For instance, what if the cascaded system had've had $(s-2)$ factors rather than $(s+2)$ factors that cancelled (ie. one of the subsystems had a pole at $s=2$). In such a case, would it not be wildly incorrect to recover the ODE above since then we are ignoring an internal characteristic mode that blows up $(e^{2t})$?


texts: I have the Fourier Transform of a signal. How can I determine the Nyquist Rate of that signal? What is the formula for that?


texts: Calculating poles/zeros from 4+ order transfer function

texts:
I need to calculate the poles and zeros from minimum 4th order (max ~20) filter coefficients. I have limitations that prohibit the companion matrix/eigen decomposition method used by roots() in Matlab.

I can see many methods for general iterative root solving, but would really appreciate some advice whether there are any methods that work well for digital filters.


texts: With respect to the dirac delta function, what does $\delta [x + 2A, y - A]$

texts:
I have come across the following function and I'm trying to figure out what it really is.  I have this function:

$f[x, y] = \delta [x + 2A, y - A] - \delta [x - A, y + A]$.  The $\delta$ is representative of the dirac delta function.  This function is used as input to some finite impulse response filter.

What does this function really mean?  Is the first part of the function 0 everywhere except $-2A, A$ and the second part 0 everywhere except $A, -A$ ?


texts: Does a matched filter give a gain in SNR?

texts:
I am trying to know if filtering a sequence with a matched filters gives a gain equal to the sequence length, he are the calculus, is it correct ?

In my model $s$ is a deterministic complex sequence (lets say RRC pulse shape samples for a QAM modulation), $v$ is a complex white circular gaussian noise.

If it is correct does it mean that there is an "inband" filter which can enhance $SNR$ ?!

\begin{align}
s &= [ s_0,\dots,s_{ N-1} ]^T\\
v &= [ v_0,\dots,v_{ N-1} ]^T \\
x &= s + v \\
R_v &= \mathbb{ E } \left[ v v ^ H \right] = \sigma ^ 2 I \\
\textrm{SNR}_1 &= \frac{ \mathbb{ E } \left[ \lvert s \rvert ^ 2   \right] } { \mathbb{ E } \left[ \lvert v \rvert ^ 2 \right] } = \frac{ \mathbb{ E } \left[ s ^ H s  \right] } { \mathbb{ E } \left[ v ^ H v \right]} = \frac{ s ^ H s } { \mathbb{ E } \left[ v ^ H v \right] }  = \frac{ s ^ H s } { \mathbb{ E } \left[ \Sigma \lvert v_i \rvert ^ 2 \right] } = \frac{ s ^ H s }{ N \sigma ^ 2 }\\
h &= s \\
y &= h ^ H x = h ^ H s + h ^ H v = s ^ H s + s ^ H v  \\
\textrm{SNR}_2 &= \frac{ \mathbb{ E } \left[ \lvert h ^ H s \rvert ^ 2   \right] } { \mathbb{ E } \left[ \lvert h ^ H  v \rvert ^ 2 \right] }
  = \frac{ h ^ H s s ^ H h }{ \mathbb{ E } \left[ h ^ H v v ^ H h  \right] }
     = \frac{ h ^ H s s ^ H h }{ h ^ H R_v h } = \frac{ s ^ H s s ^ H s }{ s ^ H \sigma ^ 2 I s } = \frac{ s ^ H s }{ \sigma ^ 2 } \\
     \frac{\textrm{SNR}_2 }{ \textrm{SNR}_1 } &= N
\end{align}


texts: Monotonic, symmetrical soft-clipping polynomial

texts:
Okay, I intend to answer this question myself, but I would like other folks to pipe in on it.  The issues or requirements are this:

Input to soft-clipper is $x(t)$, output is $y(t)$.
Output $y(t)$ limited in magnitude to $1$: $|y(t)| \le 1$
Memoryless. $y(t)$ only a function of the present $x(t)$, no past values.  So we may as well call the input $x$ and the output $y$. And we will call the memoryless soft-clipping function: $ y \triangleq f(x) $.
$f(-1)=-1$ and $f(1) = 1$.
$f(-x) = -f(x)$, an odd-symmetry (about $x=0$) function (so polarity of the input is irrelevant).
Monotonically increasing function: $f'(x) \ge 0 \quad \forall x$
For $|x| \le 1$, then $f(x)$ is a finite-order polynomial (so the frequency components generated are finite in frequency) and an odd-order polynomial (so all even-order terms are zero because of the odd-symmetry).  Let's call that odd order $2N+1$ where $N\in \mathbb{Z} \ge 0$.  Being a finite-order polynomial means a limit to the width of spurious frequency components generated which tells us what oversampling ratio suffices.
$f(x)$ to have as many derivatives equal to zero as possible where $x=\pm 1$.  So at $\pm 1$, this splices to a constant function ($f(x) = \operatorname{sgn}(x)$ for $|x| \ge 1$) with as many derivatives being continuous as possible.

So what is $f(x)$?

$$ f(x) = \begin{cases}
-1 & x \le -1 \\
\sum\limits_{n=0}^{N} a_n \, x^{2n+1} \quad & -1 \le x \le +1 \\
+1 & +1 \le x \\
\end{cases}$$

$$ $$

$$\begin{align} y(t) & \triangleq f\big(x(t)\big) \\
& = \sum\limits_{n=0}^{N} a_n \, \big(x(t)\big)^{2n+1} \\
& = x(t)\sum\limits_{n=0}^{N} a_n \, \big(x(t)\big)^{2n} \quad |x(t)| \le 1
\end{align}$$

$f(x)$ is continuous everywhere and as many derivatives as possible are continuous everywhere and the only discontinuity in any derivative is at $x = \pm 1$.

What are the odd-order polynomial coefficients $a_n$?


texts: Numericaly computing an integral in MATLAB

texts:
I have the following integral

$$\int_{-\infty}^{\infty} g(t-\tau) \, \mathrm d\tau$$

where $g(\cdot)$ is a known function, e.g., $g(t) = \sin(t)$, and I want to evaluate its value with MATLAB. How should I do this in discrete form? Thanks.


texts: Smoothed image produces contour-like ripples, yet it disappears in saved image?

texts:
I smoothed a grayscale image using a size 10 gaussian filter and display it in Octave (imshow) which resulted in ripples which appears to be contours for image intensity.

Yet if I save the image (imwrite), I get a proper smoothed image.(The ripple image above was taken as a screenshot.)


texts: How to Implement and Minimize Maximum Likelihood Expression in MATLAB

texts:
I need to check if the estimation algorithm has converged or not. I am using the Maximum Likelihood estimation method. For convergence check, we see if the log-likelihood has reached its maximum value or not. But I am having difficulty in implementing the log-likelihood expression.

As an example, I am estimating the model parameters of a Moving Average model of order d =3 expressed in Eq(1). The known coefficients are h = [1    0.45   -0.2].The pdf is given in Eq(2) and the log-likelihood in Eq(3).

I have tried the following steps. For 9 different values of the channel coefficients, I calculate the likelihood expression. I then iterate 9 times to get 9 different values. After the vector of likelihood values are obtained, I should take the log of the values of likelihood and plot. However, I am stuck in the plot as I cannot understand how to show the maximum for the coefficients of the channel which is in a vector.

Question 1:  I am stuck in implementation in the step. How do I plot the curve and show at which value of $H$ the curve reaches its maximum value?

Question 2: I want to add noise of 5dB using the awgn() function as follows

 y=awgn(y,5,'measured');

But in that case I cannot understand what value I should take for the standard deviation of the noise, s?

 clear all

N=256; %Number of Samples to collect

H_est = [1.07648398119767,0.469220009437168,-0.0245459792881367;
1.01925235521823,0.449715366555421,-0.0287160909270392;
1.05903405523260,0.568022131010140,0.0261264158909992;
1.03693902020750,0.421708336037157,-0.102062737844706;
1.07969493608754,0.481505366768731,-0.0505480628598221;
1.05948950849694,0.450984209695615,-0.110305164694739;
0.933641686721610,0.310664477120195,-0.0803322291311877;
1.13980384345926,0.530800706049278,-0.0314494051827266;
0.949926376465499,0.498110774466918,0.0473466083388137];
 %Generating source information signal
z = rand(1,N);
h = [1    0.45   -0.2]; %true known channel impulse response

y = filter(h,1,z);
y=y+randn(1,N);
s=1; %Assume standard deviation s=1
tol = 1e-4;
maxiter = 9; %because there are 9 values of the channel coefficients
llh = -inf(1,maxiter);

 llh=zeros(1,9); %Place holder for likelihoods

for iter=1:9
%Calculate Likelihoods for each parameter value in the range

 H = H_est(iter,:);

     Hz=  filter(H,1,z);
        llh(iter) = exp(-sum((y-Hz).^2)/(2*s^2));  %Neglect the constant term (1/(sqrt(2*pi)*sigma))^N as it will pull %down the likelihood value to zero for

      %  increasing value of N
        if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter-1)); break; end   % check likelihood for convergence
  end
        [maxL,index]=max(llh); %Select the parameter value with Maximum Likelihood
        display('Maximum Likelihood of H');
display(H_est(index,:));

plot(1:9,log(llh));


texts: Spectral Analysis of a Time Series with Missing Data Points

texts:
I use a PC to record time series of some physical property. The problem is that, for some reason, I did not record the time series as a whole, rather I record first segment, then second, third, etc. Each segment is 500,000 points, representing ~30 min. However, they are not strictly consecutive: between each segment and the other there will be a 'lag', or missing points, equal to a few seconds.

Since the number of missed points is too small compared to the total number of samples for each segment, can I simply concatenate the segments and treat them as a single time series? At the end I will be interested in doing FFTs, histograms, etc.


texts: Is it only me or matlab's periodogram's function is confusing?

texts:
I am really confused with matlab's periodogram.

First of all, it says in the documentation that its unit is "dB per unit frequency". What does that really mean? I just don't quite get it.

Second of all, I'm trying to evaluate an ADC. I tried to measure its linearity but that won't work because when I reduce the level of the input signal by ,let's say, -3dB I excpect to see -3dB reduction in the highest peak of the spectrum, but that's not the case.

Last but not least, with this "dB per unit frequency" when is the highest peak actually 0 dB? It should reach that at a point, right? Because as far as I know dB is always a ratio of some kind. How can I modify it as to give me the values in dBFs for example? so I could get a sense of what my data look like.

Sorry for the long post, Mathworks' Periodogram just doesn't say it clear.

Thank you all.


texts: Designing a Chirp in the Frequency Domain

texts:
I am trying to create a chirp signal in the frequency domain, by simply setting the magnitudes of each frequency bin to one, and then setting the phase delay function such that the resulting time delay increases for each frequency bin. I am trying to calculate the length of the resulting chirp from the phase delay function.
If the time delay for each frequency, $td(f)$, equals the period multiplied by the number of cycles introduced by the phase delay, and $\phi(f)$ is the phase delay in radians as a function of frequency, then
$$
td(f) = \frac{\phi(f)}{2\pi}\frac{1}{f}
$$
And if we define $\phi(f)$ as,
$$
\phi(f) = 2\pi nf^2
$$
Then it follows that the time delay, $td(f)$, is given by,
$$
td(f) = \frac{2\pi nf^2}{2\pi}\frac{1}{f} = nf
$$
And it seems obvious to me that the maximum time delay would occur when f is at its maximum $\frac{fs}{2}$, where $fs$ is the sampling frequency.
However, when I try and synthesise the signal (in Matlab),
% Create frequency axis (0 to nyquist)
fs = 44100;
fAx = 0:1:(fs/2 - 1);

% Define n so td = nf = 0.25s, and define Phases
n = 0.25/(fs/2);
phi = -2*pi*n*fAx^2;

% Calculate the maximum time delay
td = n*(fs/2);

% Calculate the phase of each frequency up to nyquist
chirp_fft = exp(1i*phases);

% Create Hermitian symmetric signal and ifft
chirp_fft = [chirp_fft(1:end-1), 0, flip(conj(chirp_fft(2:end-1)))];
chirp = ifft(chirp_fft);

% Create time axis and plot chirp
tAx = 0:1/fs:(length(chirp) - 1)/fs;
plot(tAx, chirp)

So by my calculation, the chirp should have a maximum delay of 0.25 seconds, however, when I plot it:
It is 0.5s long. When I try this for other values, the resulting chirp is consistently twice as long as my calculation says it should be. Practically this is not a problem, but I'd really like to understand where I've gone wrong! Where is my missing factor of two?


texts: Salt and Pepper impulse Denoising opencv

texts:
I am creating a generic method to work on salt and pepper noise and variants. The example images are as shown below :

I tried few methods, such as

    Median filter from scipy
    Selective Adaptive Median Filter by Jayanta Das et al.

The closest result was on Image 3, with Median filter, giving the closest result to the original image with no noise.

These are my following doubts :

Can we consider these noises as salt and pepper noise.? Is there something else that I am missing?
What could be the better suggested method? Currently, I am planning to implement Switching Median filter by Pei-Eng et al, but I would like to know if this could be the right track.

The originals that I am trying to get closer to:


texts:
I wanna experiment k-means clustering method on different kind of images, so I am trying to find different kind of images used in image segmentation such as MRI images.I want to gather some more categories.
Any suggestion would be gratefully appreciated.


texts: Extract a square-shaped signal which lies in a specific range

texts:
I need to extract (estimate) a square-shaped signal from the main signal. This square-shaped signal should be located within a specific range, in fact all we have to do is to estimate its height. A square-shaped signal is not periodical, it's length can vary from 1/12 to 2/3 of overall length.

I've read some literature about BSS and ICA but I feel like this problem should have a very simple solution.


texts: Is a wavelet-based correlation measure worth any additional computational overhead?

texts:
I have used both correlation and coherence as measures of correlation between signals.  I was thinking that a time-frequency approach would give me the best of these worlds.

My question is whether this extra data adds enough to the overall picture of the signal to justify the increased computational cost associated with doing the wavelet transforms as part of the calculation?

Reference: an ArXiv paper: "A cross-correlation technique in wavelet domain for detection of stochastic gravitational waves" by S.Klimenko, G.Mitselmakher, A.Sazonov


texts: Modulation scheme for a noisy channel with a complicated transfer function

texts:
I want to build an audio modem with a speaker in a room and a smartphone. The channel is quite bad:

Far field, so there is potentially significant reverb
Typical residential background noise
phone microphone quality, phone case, and room transfer function can be all over the place
microphone handling noise and hand partially blocking the mic

While the transfer function is unknown and can be quite variable (+/- 20 dB for narrow band features), it's safe to assume, that the channel is fairly constant during transmission.

On the positive side, it can be slow, I only need to get a handful of bits across. So I'm looking for something that's better than morse code, but doesn't require a locally flat transfer function.


texts: Why convolution in time and frequency domain are not equal?

texts:
I was writing this Jupyter Notebook. Then, I unexpectedly faced two visually similar output, but their numerical results were very different.

This is the code I used for convolution in the time domain:

filtered_k1 = sp.signal.convolve2d(img, k1, mode='same')
filtered_k2 = sp.signal.convolve2d(img, k2, mode='same')

def FT(x):
    return np.fft.fft2(x)

def iFT(x):
    return np.fft.ifft2(x)

def pad_kernel(k, x):
    x_h, x_w = x.shape
    k_h, k_w = k.shape
    padding = [[0, x_h - k_h], [0, x_w - k_w]]
    return np.pad(k, padding, mode='constant', constant_values=0)

k1_big = pad_kernel(k1, img)
k2_big = pad_kernel(k2, img)

K1_BIG = FT(k1_big)
K2_BIG = FT(k2_big)

filtered_K1 = np.real(iFT(IMG * K1_BIG))
filtered_K2 = np.real(iFT(IMG * K2_BIG))

diff1 = filtered_K1 - filtered_k1
diff2 = filtered_K2 - filtered_k2

The RSS on diff1 and diff2 was not residual numerical error, as I expected.


texts: Keypoint orientation with the oFAST from Rublee's ORB

texts:
I've just read the papers from Rublee et al. (ORB: an efﬁcient alternative to SIFT or SURF) and Rosten et al. (Machine learning for high-speed corner detection) (=fast-detector). Rublee tries to explain the oriented-FAST-Detector an its improvements.

But there is one thin I can't understand. He assumes that a detected corner is not the real centroid of this corner. He uses the image moments to determine the real centroid. He then computes the atan2 to get the orientation of the intensity patch.

The Question is: What will he do with this orientation? I don't see any sentence about a shift of a possible corner-point using the orientation and den centroid.


texts:
I am asking specifically about discrete sound signal, non dithered 16 bit 48000 Hz pcm/wav format.

My question is, how fast can you change sinewave frequency before it distort? Imagine a frequency sweep,it starts at 3000 Hz and ends at 6000 Hz, 3000 Hz full cycle sinewave at 48k sampling frequency is 16 samples long, 6000 Hz is half that, so 8 samples long...

If we do the sweep from 3k to 6K long duration,lets say one second,it will cleanly sweep from one point to another,but what is the precise limit,the threshold of how fast/short we can do this sweep before it becomes distorted?

Obviously we know its at least as long as single cycle of the lowest frequency we want to sweep,so its 16 minimum becose remember 16 samples is minimum lenght of single 3000hz cycle,but is it precisely 16 or is it longer? maybe one full cycle of starting frequency sample length(16) + same for final highest frequency(6000hz = 8) so 16 + 8?

I believe its also called "sweep rate",and sinewave frequency sweep is called "chirp",I want to know precisely the mathematics behind this to figure out where exactly this threshold is.

Also does this change depending on loudness/amplitude level? if we do the 3000 Hz to 6000 Hz sweep -6 db instead of maximum amplitude,would it make difference? considering signed 16 bit is 32768 amplitude levels,-6 db would be 16384.

I want some kind o mathematical formula that would precisely tell me exactly how many samples do I need to make shortest sweep possible at specific loudness and frequency points.So for example 3000 -> 6000 full amplitude sweep? atleast XXX samples needed! 375hz to 12000hz -12db? that be XXXXX samples minimum,you get the idea...

remember,we are talking about sweeping sinewave,not waveforms like squarewave that contain harmonics

edit
Definition of distortion,when the sinewave stops being sinewave,its simple as that.For example,if you module amplitude of sinewave slowly and only one way (up or down ) then it will always be just pure clean non-distorted sinewave,but if you start changing its amplitude too fast,so its going up and then changes direction downwards per single cycle of that sine,harmonic distortion will grow from that sine,so you will have multiple tones at same time for example 1000hz sine will get 2000hz,3000hz,4000hz... overtones/partials/harmonics all from amplitude changes

same things happens with freqency change instead of amplitude change,you do it too fast,and new sinewaves(hamonics) will appear at same moment on top of the original moduleted sine.


texts:
You can recommend the literature or tell you what the Fresnel Transform does to the image.
I realized that this describes the passage of a light wave through an optical element, but I need to understand more.


texts:
suppose I have a 4 seconds signal captured at $f_s$ = 1024Hz, which gives 4,096 data array. What is the different between these 2 FFT methods (MATLAB code)

[1] FFT the entire input with 4096-points FFT

N = length(input);
Y = fft(input)/N;

[2] Divide the input into 8 segments (4096/8=512 data point each), and feed each segment to 512-points FFT, and then average the output:

N = length(input);
Y = zeros(1,512);
for i = 1:8
  segment = input(512*(i-0)+1 : 512*i);
  Y = Y + abs(fft(segment,512))/512/8;
end

Question #1: Are these 2 methods equivalent?

Question #2: Is the frequency resolution for both methods equal to $\frac{f_s}{N}$ = 0.25Hz per bin?


texts: Complexity: Number of MAC operations on a DSP

texts:
Suppose a complex $(N \times N)$ matrix $\mathbf{T}$, multiplying a complex vector $\mathbf{x}$ with length $N$.

On a given application that I'm studying, which is irrelevant for this discussion, I can get the same results of $\mathbf{Tx}$ by the following operation:

\begin{align}
\mathbf{W_1} \Re\{ \mathbf{y} \} + \mathbf{W_2} \Im\{ \mathbf{y} \}
\end{align}
where are $\mathbf{W_1}$ and $\mathbf{W_2}$ are complex $(\frac{N}{2} \times \frac{N}{2})$ matrices, $\mathbf{y}$ is a complex vector of length $\frac{N}{2}$, $\Re$ denote the real part and $\Im$ the imaginary part.

I'm trying to asses which one requires more operations in a real implementation.

My analysis is that the first ($\mathbf{Tx}$) requires $N^2$ MACs, while the second requires

$ 2(\frac{N}{2})^2 = \frac{N^2}{2}$ MACs, plus $\frac{N}{2}$ sums.

I'm not completely sure about this, neither about how a DSP deals with complex numbers.
Am I making any mistake?

I'd appreciate your help.


texts: What is the procedure to obtain a 3D face dataset? Can anyone suggest a good 3D face dataset?

texts:
Project Description:

Input:

Frontal Face Image
Expression
Angle
Size

Details:

I have to convert the input frontal face image into 3D face and simulate the given expression (smile, cry, neutral, ..); rotate by the given angle and re-size to the given size. Again convert back into 2D face after the changes.

After googling, I found that Morphable models is a good algorithm to start with; I want a 3D face database to start & complete the implementation.

I found some 3D databases list on www.face-rec.org. Can anyone help me with buying (pricing, procedure ....) the database, and how to choose the database?

Thank you in advance..


texts: Flat fading criterion of OFDM subcarrier spacing

texts:
I have always thought that the OFDM subcarrier spacing $\Delta f$ is chosen such that

not too small because Doppler spread can destroy subcarrier orthogonality
not too large to avoid Cyclic Prefix (CP) overhead because OFDM symbol period $T_u = 1/\Delta f$ needs to be much larger than $T_{CP}$ which, to avoid OFDM symbol ISI, must be larger than multipath delay spread $\tau_m$ which depends only on the given physical environment.

The second criterion implies that $T_u = 1/\Delta f \gg T_{CP} > \tau_m \implies \Delta f \ll 1/\tau_m \sim B_c$ where $B_c$ is coherence bandwidth, this is flat fading definition. I means that flat fading is a consequence of the design criterion, not a criterion itself.

But in this wikipedia article Wikipedia's Fading Article, it is said that

  Since different frequency components of the signal are affected independently, it is highly unlikely that all parts of the signal will be simultaneously affected by a deep fade.Certain modulation schemes such as orthogonal frequency-division
  multiplexing (OFDM) and code division multiple access (CDMA) are
  well-suited to employing frequency diversity to provide robustness to
  fading. OFDM divides the wideband signal into many slowly modulated
  narrowband subcarriers, each exposed to flat fading rather than
  frequency selective fading.

It seems that flat fading is the design objective.

If a portion $B_c$-wide is in deep fade, all subcarriers in that $B_c$ are in deep fade, I am not sure that increasing diversity is the purpose. I mean : could  50 errors per 100 bits be better for channel coding than 5 error per 10 bits?

I understand that $\Delta f < B_c$ facilitates channel estimation in frequency domain by adding at least one pilot per $B_c$. Thus my question : is channel estimation the only reason to choose $\Delta f$ as small as possible ?

Now I look at the time domain, $\Delta f > B_c \implies T_u < \tau_m$. As being visualized in figure below, nothing bad happens. We can capture all the energy from all physical paths. Is there any explaination if we look only at the time domain?

Update : I am doing some math to see how the flat fading condition leads to single-tap equalizer.

Let $N$ the number of subcarrier, or FFT size (no zero subcarrier)
$T_s$ is sampling period $\implies N\Delta f = 1/T_s$ or $\Delta f T_s = 1/N$

The signal at TX after IFFT :

$$x(t) = \frac{1}{N}\sum\limits^{N-1}_{k=0} x[k]e^{j2\pi k \Delta f (t - N_{CP}T_s)}\text{ with }0 \leq t < (N_{CP} + N)$$

The received signal from $L$ multipath propagation:

$$y(nT_s) = \sum\limits^{L-1}_{l=0}h_l(n T_s)x(n T_s - l T_s)\text{ with }0 \leq n \leq N_{CP} + N-1$$

Take the part of FFT window as in the figure above, $N_{CP} \leq n \leq N_{CP} + N-1$. Set $m = n - N_{CP}$ and suppose that $h_l(t)$ is time-invariant in $0 \leq t < (N_{CP} + N)$ so that the time index of $h_l(t)$ can be dropped (this is the underspread assumption which is valid for typical channels) :

$$\begin{align}
y[m] &= y(m T_s) \\
&= \sum\limits^{L-1}_{l=0}h_l \sum\limits^{N-1}_{k=0} x[k]e^{j2\pi k \Delta f (m - l)T_s}\\
y[m] &= \sum\limits^{L-1}_{l=0}h_l \frac{1}{N}\sum\limits^{N-1}_{k=0} x[k]e^{j2\pi k (m - l) / N}
\end{align}$$

with $0 \leq m \leq N - 1$

Note that $T_s\Delta f = 1/N$.

Take FFT of $y[m]$:

$$\begin{align}
z[k_0] &= \sum\limits^{N-1}_{m=0} y[m] e^{-j2\pi k_0 m / N} \\
&= \sum\limits^{N-1}_{k=0}x[k] \sum\limits^{L-1}_{l=0} h_l e^{-j2\pi k l/N} \frac{1}{N} \sum\limits^{N-1}_{m=0} e^{j2\pi (k-k_0)m/N}
\end{align}$$

This orthogonal property assures that:
$$R(k) = \frac{1}{N} \sum\limits^{N-1}_{m=0} e^{j2\pi (k-k_0)m/N} = \delta(k_0)$$

Thus
$$z[k_0] = x[k_0] \times \sum\limits^{L-1}_{l=0} h_l e^{-j2\pi k_0 l/N} = x[k_0] \times H(k = k_0)\text{,}$$
where $H(k)$ is DFT of channel impulse reponse $(h_l, 0 \leq l \leq L-1)$ : is this the desired one-tap equalizer ?

Could someone tell me in which step I did use the condition flat fading $\Delta f < B_c \sim 1/\tau_m$ to come to single-tap equalizer model ?


texts: Deriving the Langrangian interpolation polynomials in Cook-Toom convolutions

texts:
I'm working through Blahut's 'Fast Algorithms for Signal Processing'. Trying to develop an intuition for the Cook-Toom algorithm for convolutions as used by Lavin and Gray in their Winograd paper for convnets.

At page 151 (chapter 5.2), Blahut writes that the following is the final result of the Cook-Toom algorithm:

$s(x)=s(\beta_0)L_0(x)+s(\beta_1)L_1(x)+s(\beta_2)L_2(x)$

I'm a bit confused about why one has an $s(x)$ and then has an $s$ instead a few lines later like so:

$s=C[Bg·Ad]$

But that is secondary. I'm guessing that's to do with us no longer treating $s$ as a polynomial after the inverse transform with $C$? But I though convolution was equivalent to polynomial multiplication. Still not entirely clear on that.

My main question is, how are the following Lagrange polynomials derived? For reference, the interpolation points are $\beta_0 = 0, \beta_1 = 1, \beta_2 = -1$:

$L_0 = -x^2 + 1,$

$L_1 = x^2 + x,$

$L_2 = x^2 - x.$

These are stated, but don't follow obviously from the rest of the text. My understanding of Lagrange polynomials is somewhat limited. Any clarification would be very much appreciated!


texts: Regarding phase and frequency components in the signal

texts:
As I am not a signal processing student, I have limited understanding of the concepts.

The general formula for sine wave is

$$s(t)= A \sin(\omega t-p)$$

where $A$ is amplitude $\omega$ is frequency and $p$ is phase.

Is there any possibility that the phase of the signal remains constant even if the frequency of the signal is changed?
If we observe a real-time signal the signal can be a periodic sine wave signal but not smooth curve as it is accompanied by noise. In such cases how the phase of the signal is varied? Can the noise effect the phase of signal?


texts: Is it possible to use a Volterra series to generate subharmonics?

texts:
Taylor series can create harmonics of the frequencies in an input signal. I'm wondering if it is likewise possible to use Volterra series to create sub-harmonics of the frequencies in that signal.

Some degree of finesse is required in posing the question:

Taylor series don't only create harmonics; there's intermodulation distortion as well. Likewise, it's ok with me if a suitable Volterra series doesn't only create sub-harmonics, but also creates something similar to intermodulation distortion, or even some degree of harmonic distortion.

However, I would like to exclude trivial cases where the input signal gets so wrecked that it ends up looking like white noise, which would obviously create sub-harmonics (along with every other frequency you can imagine!). I'm not sure how to formalize this requirement precisely, other than to say that I hope it's clear what I'm driving at.

Here's one way to formalize the behaviour I want:

Given an input $\cos(\omega t)$, yields an output containing $\cos(\frac{\omega}{n} t)$ for
some fixed $n$, ideally a natural number
Given a general input, a sub-harmonic is generated corresponding to every frequency in the input
Other "intermodulation distortion"-type artifacts are allowed, "within reason"

The answer can be expressed for either discrete or continuous signals.


texts: Using voice recognition algorithms to reduce noise in financial time-series?

texts:
Is it possible to use voice recognition algorithms to reduce/remove the noise in financial time-series?
Does such an approach even make sense?

If "yes", how could that be done technically?


texts: How does time shift correspond to phase change in a discrete signal?

texts:
I was watching this video where the presenter remarks:

  For a discrete signal, time shift corresponds to phase change in a discrete signal but not vice versa.

I was trying to figure out how this could be proved. For example, if our discrete signal is

$f[n] = A \cos(\Omega_o n + n_o)$ where $ n,n_o \in Z$

If the signal undergoes a time shift $n'$, we have

$ A \cos(\Omega_o (n + n') + n_o) = A \cos(\Omega_o n + \Omega_o n' + n_o)$

This would correspond to a phase shift only if $\Omega_o n'$ is an integer. But that means $\Omega_o$ needs to be an integer, which is not always the case.

What is the correct explanation here ?


texts: Baseline correction for continuous EEG signals

texts:
I'm new to the topic so forgive the trivial question.

I am trying to understand data preprocessing for EEG signals. As I understand, after recording continuous EEG signals, they should be band-pass filtered and segmented into epochs (so the window of time that is relevant to the experiment). After that, these epochs should be baseline corrected. Can anyone explain what this concept of baseline correction means?


texts: Average identical sweeps or lowpass filter sweeps first?

texts:
There are 5-10 replicate sweeps of measurements that are meant to be averaged. I also typically lowpass filter the signal at 2000Hz. The target signals are in the 0-100Hz range, and have, for lack of a better term, varying variability (the 5-10 sweeps are usually very similar to each other, but sometimes they are pretty different, but averaging out the variability is fine since it's random).
My question is whether I should (1) lowpass filter the 5-10 sweeps individually before averaging them, or (2) average the 5-10 sweeps and then lowpass filter the resulting mean.
Is there a general rule that applies to this? I am new to signal processing, but I'm aware that order of operations matters.


texts:
I have a signal with bandwidth between 0-2KHz.
If i would like to know the signal's frequency spectrum energy, then I use Parseval's theorem and calculate it over time.
Now it's my question:
What if I would like to know only part of the spectrum energy, let's say 1.5Khz-2Khz, Is there a low cost way to do it, instead of FFT the signal and then sum the relevant frequency band?


texts: Removing 'high-frequency' components from an array?

texts:
I am new to the stream of digital signal processing, and hence my doubt may seem extremely stupid, but:
In tools such as matlab, signals are stored in a single array. If i convolute it with an LPF array, i get another array/signal, that has certain elements attenuated, can i find out which elements these are and completely remove them from the new array.

Does my question even make sense? :/


texts:
I am currently trying to detect heads in a sequence of real-footage images and am using HOG feature descriptor and SVM as classifier.

Currently I am using Dalal's HOG implementation code in MATLAB found in this link:
http://www.mathworks.com/matlabcentral/fileexchange/46408-histogram-of-oriented-gradients--hog--code-using-matlab

Currently I am using libSVM MATLAB version found in this link:
http://www.csie.ntu.edu.tw/~cjlin/libsvm/

I prepared 350 positive training images and 1243 negative training images.

The hog feature vectors extracted from the training images are converted to libsvmFormat and inputted in the libsvm training method, to obtain a model. The hog vector length is that of 1764.

Regarding the libSVM I chose these as parameters:

-s 0  (i.e. C-SVC)
-c 1  (i.e. cost = 1)
-t 2  (i.e. kernel = RBF)
-g 3  (i.e. gamma = 3 (this is for kernel))

Regarding the HOG version, I left the cell, bin & block parameters as they were in the implementation shown in the link above.

I am using a scanning window of size 128x128 and 256x256 to scan through the whole image to detect possible heads.  At each window, the hog feature vector is extracted for each image and is inputted in libsvm predict, to test whether it should be classified as a head or not.

However, after doing all the above, I have a numerous amount of false negatives and can't figure out what I am doing wrong.

Can someone experience please offer some advice on what is possibly wrong?  I really need to figure this out please.  Much appreciated!


texts: Spectral Property of Bandlimited Quantized Signal

texts:
recently I meet a problem.
I generate a random Gaussian signal $x$ with mean $0$ in MATLAB. And use a $4$th lowpass Butterworth filter to acquire $x'$. Then I use the binary rounding quantization method to get the quantized signal $x''$. And I find that the Fourier spectrum (FFT) of $x''$ contains some high-frequency information. That means some energy leaks to high-frequency regions after quantization.
Also, I find the energy in high-frequency areas is related to the bandwidth of the filtered signal $x'$. When the bandwidth of the filtered signal $x'$ is relatively low,  the energy of the high-frequency areas of the quantized signal $x''$ is also lower. I find that some historical references about the assumption of the whiteness for quantization noise. But it seems that this is not the case.
Do you all meet this question? Could anyone give me references for this problem? Thanks in advance!


texts:
I have many close-range images of a single wind turbine blade (taken from different positions, same angle):

I need to separate the blade image from background like this:

The problem arises in image alignment application where the blade image moves about large amount while the background does not (parallax effect). The backgroud also consists of more pixels so it is prefered for alignment while the blade is discarded.

Another problem is that the blade surface is very smooth and there are few to none interest points to detect.

Is it for example possible to detect one or more curves to isolate boundary of the blades and then possibly align the curves instead of interest points?

Another approach would be to use graph cut to extract the precise edge and convert it to spline curve - however, how to align such curves? Maybe 3D surfaces based on intensity data can be aligned rather than 2D curves.

Any more ideas?


texts: How can I convert time-value pairs to a signal with a constant sampling frequency?

texts:
I have a time series consisting of [time, value] pairs. The series is ordered by time but the difference between two consequent times is not constant. For example:

time, value
[0.1,  5.2]
[0.23, 6.9]
[0.3,  7.2]
[0.31, 4.8]
...

What are the methods to convert such time series to a vector with a constant sampling frequency? This will probably involve some interpolation of values that were nonexistent before, but that is ok for me.

Is there any matlab function for doing so?

The final goal is to be able to apply filters and FFT on the original timeseries.


texts:
I want to design a Filter with a variable Cutoff Frequency and Resonance. However in the source I'm using there is a Bi-Quad Filter with two IIR and two FIR Filters designed like this :

Furthermore I found that I can compute the coefficients for the different Filter types as stated in this table :

Is there a way to alter those coefficients that i can compute the samples with a variant FIlter Resonance?

EDIT Since I'm designing this filter for a software synthesiser, with resonance I mean the typical boost or reduction of the cutoff-frequency in a variable manner.


texts:
I want to run a numerical simulation of a system (eg first order low pass filter) that is stimulated with a photocurrent-like signal. My preferred language is Python.

I'm failing to understand how to generate a proper time series for the photocurrent with a consistent noise power, especially when the photocurrent becomes very small.

Let's assume:

simulation interval time is dT
the mean photocurrent is I
q is the elementary charge
T is the time series length

In my understanding there are two main ways to model the photocurrent time series.

A) I calculate the average number of electrons and use this number to generate a poisson time series that I convert back into current

N = I * dT / q # mean number of electrons in interval dT
N_ts = np.random.poisson(N,int(T/dT)) # electrons time series
I_ts = N_ts * q / dT # resulting photocurrent time series

The problem I have with this is that if I becomes small enough then average number of electrons per dT is very low and when I draw samples from the distribution it behaves like a quantum current (some samples will have 0,1,2,3,.. electrons). Is this quantization a corrected model or is it a limitation of it? What is the meaning of dT in this case? Is it some sort of averaging? Should the dT be increased to avoid the quantization? If I increase dT, how much should it be increased? Does a change in dT change the noise power of a Poisson time series? How is dT related to the noise bandwidth of the signal? How can I rescale the noise after a dT change to maintain the same input noise power?

B) A second approach could be to avoid the photon counting and directly model a photocurrent time series drawing samples from a normal distribution with $\mu_I = I$ and $ \sigma_I = \sqrt{2qIB} $
However, what is the correct bandwidth to use in this case? Is it 1/(2*dT) because I'm assuming the sampling frequency is two times the noise bandwidth? A problem with modeling using a normal distribution is that if dT=1E-6 and I<1.602E-13 then $\sigma_I=\mu_I$, which means that if I draw samples from such distribution I might end up drawing negative currents. The normal distribution is not a good approximation of a Poissonian distribution for small currents. Should I use a lognormal distribution instead?

Are these approaches fundamentally correct? How do they differ from each other? Is there an equivalence between the two? Does Physics come actually into play with quantization when considering very small currents or is it an artifact of the modeling?


texts: What Techniques Are Used for Signal / Image Denoising (Example Image Provided)?

texts:
I have an image that looks like this (it might appear low res in a browser because it is 16 MB). This image was taken with a scanning electron microscope, but because the equipment is outdated there is a lot of noise within the image.

What I would like to do is develop a program (in Java) that would take this image as an input and try to remove the noise as much as possible to produce a high quality image.

However, being a beginner in the field of image/signal processing, I have no idea where to begin researching techniques. I was wondering if you guys had any good starting points for me to look into, or perhaps de-noising algorithms/techniques that I could research to potentially help me reach my goal.


texts: Cross power spectral density (CPSD) comparison between signals

texts:
I want to compare various signals with Matlab's CPSD to identify shared frequency components. (If there's a better method, I would like to know!)

I have 3 signals, A, B, and C. B is pure noise while A and C share nearby frequencies.

I noticed that CPSD's output is scaled by the power input of the signal - so for comparison between A/B/C, they should have normalized power so the CPSD amplitudes are comparable to each other (hopefully?).

The issue is that if I normalize the power, the B noise signal is boosted above the true 4 Hz signal and it looks like there's minor 5 Hz stuff for the B noise signal as well - which are both false.

How can I compare A/B/C effectively and see the relations but avoid boosted noise amplitudes in CPSD output?

I've also included a periodogram which shows the stronger 5 Hz signal but not the 4 Hz signal, for reference. mscohere gave gibberish as well. Finally, I want to keep this like "exploratory" and pretend I don't know there's a low-level 4 Hz signal that I can filter for specifically.

Any suggestions would be greatly appreciated, as I could be going in a very wrong direction or doing something incorrectly and have no idea!

t = (0:6:1080)'; %sec, some time steps
rng(1863)
A = 2*sin( 1/5*t ) + 0.25*sin( 1/4*t) + 2*(rand(length(t),1)-0.5); %some data with noise  + 2*(np.random.uniform(size=t.size)-0.5)
B = 1*(rand(length(t),1)-0.5); %some data with noise
C = 0.75*sin( 1/4.8*t ) + 0.25*sin( 1/4.05*t) + 2*(rand(length(t),1)-0.5); %some data with noise

window = hamming(110); %some window
windowPad = padarray(window,size(t,1)-size(window,1),0,'post'); %padded to t size for periodogram
Fs = 1/(1); %1/sec, time delta in freq form

pwr_A = sqrt(1/length(A)*sum(A.^2)); %estimate power of signal
pwr_B = sqrt(1/length(B)*sum(B.^2)); %estimate power of signal
pwr_C = sqrt(1/length(C)*sum(C.^2)); %estimate power of signal

[Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
Axy_AvB = angle(Cxy_AvB)*180/pi;
Pxy_AvB = abs(Cxy_AvB);

[Cxy_AvC,freqs_AvC] = cpsd(A,C,window,100,512,Fs);
Axy_AvC = angle(Cxy_AvC)*180/pi;
Pxy_AvC = abs(Cxy_AvC);

[Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
Axy_BvC = angle(Cxy_BvC)*180/pi;
Pxy_BvC = abs(Cxy_BvC);

figure(1);
subplot(2,1,1);
plot(1./freqs_AvB,Pxy_AvB);
hold on;
plot(1./freqs_AvC,Pxy_AvC);
plot(1./freqs_BvC,Pxy_BvC);
xlim([0 10]);
xlabel('Periods (sec)');
ylabel('Arb. Power');
title('Cross power spectral density matlab');
legend('AvB(noise)','AvC','B(noise)vC', 'Location','northwest');

A = 1/pwr_A*A; %normalize to a power of 1
B = 1/pwr_B*B; %normalize to a power of 1
C = 1/pwr_C*C; %normalize to a power of 1

[Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
Axy_AvB = angle(Cxy_AvB)*180/pi;
Pxy_AvB = abs(Cxy_AvB);

[Cxy_AvC,freqs_AvC] = cpsd(A,C,window,100,512,Fs);
Axy_AvC = angle(Cxy_AvC)*180/pi;
Pxy_AvC = abs(Cxy_AvC);

[Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
Axy_BvC = angle(Cxy_BvC)*180/pi;
Pxy_BvC = abs(Cxy_BvC);

subplot(2,1,2);
plot(1./freqs_AvB,Pxy_AvB);
hold on;
plot(1./freqs_AvC,Pxy_AvC);
plot(1./freqs_BvC,Pxy_BvC);
xlim([0 10]);
xlabel('Periods (sec)');
ylabel('Arb. Power');
title('Normalized Power Cross power spectral density matlab');
legend('AvB(noise)','AvC','B(noise)vC', 'Location','northwest');

[pxx,f] = periodogram(A,windowPad,512,Fs);
pxx = abs(pxx);

figure(2)
plot(1./f,pxx)
hold on;

[pxx,f] = periodogram(B,windowPad,512,Fs);
pxx = abs(pxx);

plot(1./f,pxx)

[pxx,f] = periodogram(C,windowPad,512,Fs);
pxx = abs(pxx);

plot(1./f,pxx)
xlabel('Period (sec)');
ylabel('Normalized Power');
xlim([0 10]);
title('Normalized DFT of A,B,C');
legend('A','B','C', 'Location','northwest');


texts: Impulse Invariant method for digital filter design

texts:
One of the known methods for discretizing analog filters is impulse response invariant. We get the impulse response in time domain, discretize it and then get the Z transform.
What I am trying to understand is why the freq response of the resulting digital filter has a freq response magnitude scaled by (1/T) T:sampling time?
Matlab, using c2d command, modifies it by multiplying by T so that the freq response is similar to the analog filter, but this is not the result of the Z transform I described earlier.


texts: 3D (time, scale, amplitude) plot in Continuous Wavelet Transform

texts:
I will be extremely grateful if someone could please answer this basic question.
How can one plot a 3D  (translation, scale, amplitude) plot from the Continuous wavelet transform (CWT) coefficients?
The CWT coefficient is a M X N matrix. Which of the axis corresponds to the translation, scale and Magnitude.
Thanks for your help.


texts: Comparing Welch and FFT power spectrum density

texts:
Based on previous answers from the forum, I implemented a function to compute the Power Spectrum of a 1D time series.

def pow_spect(x, fs):
    nt = len(x)
    power = np.abs(np.fft.rfft(x, nt))
    ff = np.fft.rfftfreq(nt, 1. / fs)
    power = (power ** 2) / (nt ** 2)
    power *= 2.
    power[:, 0] /= 2.
    if nt % 2 == 0:
        power[:, -1] /= 2.
    return ff, power

I compare the output of this implementation to the estimation of the power spectrum density given by scipy.welch :

nt = len(x)
ff_welch, pxx = signal.welch(x, fs, window=signal.get_window('boxcar', nt), scaling='spectrum')

Eventually, I compute ff, p = pow_spect(x, fs) and I compare : p_db_1 = 10. * np.log10(p) and p_db_2 = 10. * np.log10(pxx) for a given time series x. I noticed that all the coefficients of p_db_1 and p_db_2 are very close (up 10^(-8)) except for the first one, p_db_1[0] and p_db_2[0]. Did I make a mistake somewhere ? Thank you for your help.


texts: Is it possible to approximate given skewness and kurtosis values using an IFFT?

texts:
I have time domain data for a signal that behaves randomly but has non-random (non-Gaussian?) skewness and kurtosis values of ~0.9 and ~7.4, respectively.

The FFT of the signal shows that it's not broadband random -- it's generally flat below 100 Hz and tapers off at about -3 dB/oct above 100 Hz up to a Nyquist of 400 Hz.

I would like to create a synthesized signal that is statistically equivalent (or approximately so) to my reference signal, with nearly identical

Frequency content
RMS
Skewness, and
Kurtosis.

Additionally, I have extrapolated PSD amplitudes up to 10 kHz that I would like to add into the mix and invert to the time domain.

I have done this in the past satisfying only (1) and (2) by using, for example, MATLAB's randn() function to create random white noise (namely to create random phase angles), applying an FFT to the white noise, scaling its frequency components to match a desired FFT spectrum, and performing an IFFT to convert to time domain.  The rand() function approaches a skewness and kurtosis of 0.0 and 3.0, respectively, so my initial thought was to create white noise with a specified skewness and kurtosis.

I have tried this same procedure using MATLAB's pearsnd() function, which allows one to specify skewness and kurtosis values, but the above procedure doesn't seem to preserve the skewness and kurtosis of the scaled time domain signal.  I believe the IFFT assumes a Gaussian distribution.  Is there a work-around for this?

Thanks!


texts: What is the difference between `image normalization`, `contrast stretching` and `histogram equalization`?

texts: I believe that they all do the same task of distributing the frequency of pixel intensities to a bigger range, but can these terms be used interchangeably?


texts: How to save Fourier series approximated signal to a WAV file

texts:
I changed this Matlab/Octave code to approximate square wave by using combination of Fourier series and Fejér taper:
    % needed for Octave ------------------------------
    pkg load signal
    % ------------------------------------------------

    % Stage 1 ----------------------------------------
    figure
    N=20;                               % harmonics
    fs = 44100;
    fc = 440;
    w = 2 * pi * fc;
    r_start=0;
    r_end=10;
    r_step=1/fs;
    t = r_start:r_step:r_end;           % range
    w_sqr = square(w*t);          % Square wave [-1:1]

    plot(t,w_sqr);
    axis([r_start r_end -1.2 1.2]);
    grid on;
    hold on;
    audiowrite("squarewave.vaw", w_sqr, fs);

    % Approximation: Combined
    % Fourier series-Fejér taper (Tukey window?) -----------
    %i=1;
    %sum=0;
    %r=t;
    %for t=r
    %  for n=1:N
    %    %sum = sum + (2*sin(n*t)+sin(pi*n-n*t)-sin(n*t+pi*n))/(2*pi*n) % Fourier series
    %    %sum = sum + ((1-n/(N+1))*(-1)*sin(n*t)*(cos(n*pi)-1))/(n*pi); % Fejér taper
    %    sum = sum + ((n-2*N-2)*(cos(pi*n)-1)*sin(n*w*t))/(2*pi*n*(N + 1)); % (Fourier+Fejér)/2
    %  end
    %   Ffej(i)= sign(sum)*(1/2) + sum;
    %  i=i+1;
    %  sum=0;
    %end;
    %Ffej=Ffej';

    sum_ = sum((((1:N)-2*N-2).*(cos(pi*(1:N))-1).*sin((1:N).*w.*(0:1/fs:1)'))./(2*pi*(1:N)*(N + 1)), 2);
    Ffej = 2*sum_;
    Ffej = Ffej';

    plot(t,Ffej);
    axis([r_start r_end -1.2 1.2]);
    legend('Square', 'Fourier series-Fejér taper');

    audiowrite("approximated_squarewave.wav, Ffej, fs);

which results this plot:

I would like to save these signals to .wav files by using Matlab/Octave audiowrite() command.
What changes are needed to be done in code to save an approximation of 440Hz square wave sound sampled at 44100 Hz? My thought is that range t must be changed to support sample rate range and frequency but, how is the situation with calculation of Ffej ... or is data of Ffej directly usable in audiowrite() command?


texts: How to make series circle shift invariant

texts:
I have a series of numbers. They may be integer, floating, or binary. Length of series is constant. I need criteria which is circle shift invariant. For more information assume series S1, S2, S3, S4 and S5:

S1 = {A, B, C, D, E, F}
S2 = {F, A, B, C, D, E}
S3 = {D, E, F, A, B, C}
S4 = {A, C, D, F, E, B}
S5 = {A, B, C, D, E, G}

where A, B, C, E, F, and G are numbers. S2 and S3 is resulted by shift circle of S1. So the criteria should return same values for these three series. S4 have same number of S1, S2, and S3 but order of numbers are not the same. S5 and S1 are similar so I want to criteria of S1 and S5 be similar.

I want to use them as descriptors in image processing. In other words, I want to make a feature (like LBP, or HOG) rotation invariant. For example mean and variance can be used but I need more. Are there any suggestion?


texts: Making longer a time domain signal by adding values on its frequency domain. What am I doing wrong?

texts:
Background

Here's the thing: using software for Finite Element Acoustic Simulation I got the dataset of frequency response from a room; software works by solving the wave equation in the interval $`[f_i,f_f]`$ so results come directly in the frequency domain and we save them into a .csv file:

For decay time analysis, we're on the way to transform the data to time domain (using IFFT) and it has some tricks that you may know but I'll mention briefly:

(From Matlab's FFT documentation)

Space between 0 [Hz] and f_i [Hz] must be filled in (keeping the distance between points as in the results) and assigned zero amplitude, to let IFFT algorith work correctly
Do an amplitude "de-normalization"
"Reflect" the frequency spectra using the last frequency element as pivot, omitting the first one and conjugating this copied values
Apply IFFT, and done.

So far, we've made and proved a code that does this steps using the .csv file as entry. The test was realized with a .csv created from a .wav file and it resulted ok:

(Subplot 1: original. Subplot 2: after process)

The problem

Our analysis was made from 0.1 to 200 Hz with steps of 0.1 Hz. When we transformed the results to time domain, the signal was 10 secs long but we expected it to be at least 30 secs (since was a model similar to last's plot one). We tried to reduce the distance between points by interpolating it, as $`\Delta f = 1/T`$ say (with $`T`$ as signal duration).
However, when applying again IFFT the plot was exactly the same!

So, how can I process my data to get a longer signal in the time domain?

Hope I have been clear enough. Thanks in advance.

EDIT: showing correctly the list before first image


texts:
Is there a simple way to generate a skew-normal shape? Shouldn't be perfect, I just need a rough approximation of it to test peak-detection algorithms.
I tried generating the signal according to the continuous formula:
y = 2 * gaussian(x) * cumulativeGaussian(skewness * x)

But the problem is (I think) - we can't implement infinity. So if I pass left skewness - this won't let cumulativeGaussian to ever reach large values - it will stop somewhere in the middle of the cumulative signal.
PS: I found this answer/this example, but they talk about distributions of random variables which we can sample from. But I simply need to generate the shape with N points, I don't need a RNG.


texts: Interpretation of steady state sensor signal

texts:
I try to understand some measurement results of an accelerometer MEMS sensor. I used a mikrocontroller which is sampling an accelerometer's analog signal with the built in analog-digital-converter. Sampling rate (timer controlled task, starting the ADC) is 760 Hz, the signal is decimated to 100 Hz (cutoff @ 50 Hz) and then lowpassed to 3 Hz.

My inital goal was to have a look at the noise performance, i.e. how many ADC LSB noise is present. Therefore the µC logged the values for ~16 h and sent the data to a computer via a very slow serial interface. Transfer rate was only 1/5 Hz (otherwise a long time logging would have been complicated).

In the pictures you can see the results of three adc channels. To get an idea of the noise the first sample has been subtracted from all channels in the first plot.
There is noise -2 to 1 lsb digit in all channels (x,y,z).

My software is showing a power spectrum automatically (using an hanning window), these are the results (first plot for original data, second plot for data with 1st sample subtracted)

My question is if I can interprete something from the periodic peaks in the original signal which have a distance of 0.0096 Hz which is about 104 seconds in time domain.
And why there are no peaks in the "pure" noise signal (plot 2)? Are they simply an artifact?

Are there any deterministic components in this signal (which I might be able to correct) or is this noise stochastic?

Thank you very much!


texts: Quadratic Programming with Linear Equality Constraints

texts:
I need to solve an equality constrained minimization problem as give below
$$\min_{\textbf{w}} \mathbf{w}^TR\mathbf{w} $$
such that
$$X\mathbf{w} = \mathbf{1}$$
where $R\in \mathbb{R}^{n\times n}$ is covariance matrix (hence positive semi-definite)

$X\in \mathbb{R}^{m\times n}$ matrix with $m\gg n$,

$\mathbf{1}$ is column vector of size $m$ with all $1's$

$\mathbf{w}\in \mathbb{R}^n$ is an unknown vector.

Is there any closed form solution for this? If yes, can anyone provide it?

If no can we solve it by GD or any other algorithm?

Thanks in advance.


texts:
In my research work i have the equation y(n) = hx(n)+w(n). My understanding about this equation from MATLAB point of view are:

h has to be a single number. (for e.g., rayleigh generated).
x(n) , w(n) will be vector because they have n associated with them. So is my understanding correct or is their any other concept. Thanks.


texts: Can ICA be applied to solve BSS problem in frequency domain?

texts: I'm recently doing research on EMI, and I want to seperate the mixing signals. However, the data given is in frequency domain. As I have used ICA to solve blind source separation problem in time domain, can this be used to solve BSS in frequency domain? And how? I can not imagine how this would work. I would appreciate if you could provide me some extra materials on this subjects!


texts: What is meant by "correlation" when referring to spectral coherence

texts:
I've been reading that coherence measures the correlation between two waves as a function of frequency. I also read that difference in phase does not mean less coherence at a given frequency, and that amplitude is also not a factor. Some explanations say that its is a measure of the shift in phase. As I understand it, a shift in phase requires a shift in frequency, but coherence is the measure of correlation at a given frequency. So I guess what I want to know is this:

If I've taken the Fourier transform of two waves, and am looking to see how they are correlated in the frequency domain, what am I actually measuring?

EDIT: I'm a biologist, and my knowledge of digital signal processing is pretty basic. I'm trying to pick this stuff up on the fly. So intuitive explanations, if possible, would be greatly appreciated.


texts:
I have a question about invariances in neural networks.

In general, neural nets with enough layers can learn arbitrarily complicated nonlinear functions. Therefore, it's not hard to understand how a neural net can learn a Fourier transform or convolution: (3rd to last paragraph of this page: http://www.dspguide.com/ch26/3.htm).
Fourier transforms have a useful translation (shift) invariance property, i.e. the magnitude of the FT of a function/image remains the same after you translate it.
Other transforms also have useful invariance properties, e.g. the Mellin transform is scale invariant, i.e. the magnitude of the Mellin transform of a function/image remains the same after you scale it.

So my question is about why neural networks don't already have built in invariances. Is it purely a matter of number of neurons? If I have an enormously large number of neurons and/or layers, do I start to get built in invariances?

(E.g. as described in above link, you need 2 layers of N nodes to learn convolution of an N sample signal. Also, not every single neuron is going to be weighted in such a way to participate in learning convolution, so you'll need many more neurons than just 2N.)


texts:
I have a predicted mel log spectrogram $S$ and I want to add a noise signal $\eta$ to it.

S = log(mel-filterbank(abs/real part of(stft(signal)))

The naive way would be to approximate the linear spectrogram, add the noise signal and then apply the log mel transform again. However, it seems to me that this approximation adds unwanted artifacts which would not show up in a natural signal.

Is there a way to combine the signal without doing the approximation, given that I want the combined signal to be a mel log spectrogram too?


texts:
I have to calculate the bitrate of PAL format $\ (576*720) $. I have the following data:

$\ 25 $ frames per second
Interlaced
Interior representation is: $\ 4:2:0 $
$\ 8 $ bpp

I have tried the following formula but I am really not sure at all if it is correct because I have not clearly understood the difference between interlaced and non-interlaced video and how this difference changes the calculation fo the bitrate.

This is what I tried:

$\ (576*720)*(1+1/4+1/4)*8*25 $

This part: $\ (1+1/4+1/4) $ is due to the specific representation $\ 4:2:0$


texts:
I receive complex I and Q data from an USRP. I want to do a Hamming window LP filter on the data. However, I don't get the grasp of how I use the complex data. It is unusual because most of the time you sample in time domain and for that situation I get that you need to do a convolution with the coefficient of your filter. Should I convert the complex data in time data? or should I do a multiplication in the frequency domain?

The second part of my question concern the analyze of that complex data stream. If I want to see the frequency content of my signal, can i just take the real part and imag part of the signal to get my FFT or do I need to do a FFT over the real and imag part or the norm?


texts: How to compute the main lobe width of generalized Hamming windows (i.e. Hann, Hamming, etc.)?

texts:
Is there a way to compute the main lobe width of windows of the generalized Hamming window family (i.e. Hann, Hamming, etc.)? By main lobe I mean the first zeros left and right of the center of the Fourier transform of the windows.

I could not find any derivation in any of the standard textbooks. In Oppenheim/Schafer 'Discrete-time signal processing' (2nd Ed., p. 471), only the 'approximate width' of the main lobe is provided without proof.

Just to give you some background: Windows of the generalized Hamming window family are given by the following equation:
$w_H(n) = w_R(n) \cdot (\alpha + 2\beta \cdot cos(\Omega_M n))$ (see e.g. here), where $w_R(n)$ is the rectangular window of length $M$, and $\Omega_M = \frac{2\pi}{M}$.

The Fourier transform (FT) of these windows is given by:
$W_H(\omega) = \alpha W_R(\omega) + \beta W_R(\omega - \Omega_M) + \beta W_R(\omega + \Omega_M)$.
In this equation $W_R(\omega)$ denotes the Fourier transform of the rectangular window.

$W_R(\omega)$ is simply a sinc-function, so $W_H(\omega)$ consists of the sum of three sinc functions with centers at different positions of the frequency axis.


texts: Simple beat detection algorithm for microcontroller

texts:
For the purpose of controlling DIY disco lighting, I'm looking at implementing beat detection in a microcontroller with limited computing power (i.e. Arduino Uno). While detecting all beats would be a plus (for one, it would allow measuring BPM, which could be used to adjust the light to the style of music), it's not absolutely needed for this application since light setting must be typically changed at a much lower frequency than the beat. So what I'm really looking at is to detect "conspicuous" events in the music (beat, resume after blank, etc.) to synchronise light change, but without 100% reliability.

The constraints are the following:

Limited computing power is available. The Uno is based on a ATmega328 8-bit 16MHz MCU. If this is too constraining, I could consider moving to a Due, which is based on a AT91SAM3X8E, 32-bit 84MHz Cortex-M3 MCU.
Noisy sound input. Ideally, the sound would be captured with a small, cheap microphone. I'm not sure yet what sampling rate can be achieved with the Uno. In any case, I expect the signal to be fairly noisy. If that proved the be required, I could possibly consider implementing a proper sound input (e.g. with the rugged audio shield) but I'd rather avoid that if possible because of the constraint of hooking the music playing devices (there may be several) to the light control electronics.
Ideally, the algorithm would have a parameter to adjust sensitivity to beat event detection. That way, the sensitivity could be automatically increased to make sure events are detected even in bad conditions (low music volume, noise from conversations, etc.)

My questions are thus:

What (class of) algorithm would fit appropriately this application?
Is there any implementation available that could be readily integrated in a Arduino project?


texts: High modulation index PSK - carrier recovery

texts: How to recover the carrier of a PSK signal that has high modulation index (eg. 750 symbols per second at carrier frequency 1500 Hz)?  I guess that a PLL will not be able to lock onto it with such frequent phase transitions. Maybe upconverting or frequency multiplication of the signal before the PLL is the solution?


texts: What are the typical causes of clicks and pops in an audio file?

texts: If a device generates clicks and pops at start and end of audio file playback, even when those audio files are completely zeroed out and dc-offset corrected, what could've caused them?


texts: 2D Haar Transform Computational Complexity

texts:
I want to learn the computational complexity of the 2D Discrete Haar (Wavelet) Transform (DHWT). The number of operations (divisions, summations etc.) is the main focus. I think there are many different ideas on that. There are even papers claiming that we do not need any multiplications for it. I am confused. Any help? Thanks.

Edit: I should have added the orthogonality requirement for the Haar transform.


texts: Magnitude response and DFT normalization

texts:
Suppose I have an FIR denoted h that represents the impulse response of a system. Using MATLAB syntax (for convenience and brevity),

What does abs(fft(h)) represent?

If the answer is the system's magnitude response, then please proceed to the next question. If not, please correct me.

Is it correct to say that the value of the magnitude response of a system at a given frequency represents the gain applied to an input sine tone at that frequency?

If yes, then please proceed to the next question. If not, please correct me.

If I upsample h by, say, a factor of 2, then the values of abs(fft(h)) are doubled compared to the same computation done at the original sampling rate, at least up to the Nyquist frequency of the original rate. This seems to contradict the answer to question 2, because even though the impulse response has twice the number of samples, it represents the same system and therefore, if abs(fft(h)) represents the magnitude response (with magnitude response defined as mentioned above), then there shouldn't be a discrepancy.

These questions should reveal gaps in my knowledge regarding the meaning of a system's magnitude response, and normalization when computing a DFT (see the answers to this related question). I'd appreciate some clarity on these points. Thanks!

Update 1: Below is a figure showing an example.

The chosen example FIR is a low-pass filter. Blue is the original sampling rate (1024 samples/second). Red is upsampled by a factor of 2. The system magnitude response is a result of doing abs(fft(.)) to each of the IRs shown on the top left. The two plots below (bottom left should actually read "spectrum" and not "response") shows what the system does to an input, unit amplitude sinusoid of frequency 50 Hz.

Update 2: The confusion has been resolved in the comments section of the accepted answer. The source of confusion was the assumption that upsampling a discrete-time impulse response (as shown above) is exactly the same as sampling an analog IR at the higher sampling rate. This is not the case.


texts:
I'm doing research on image and audio watermarking, and so far I've read some literature in order to get basic introduction and good initial understanding of this topic. Right now I'm reading book "Digital Watermarking and Steganography" by Cox, et al.

My basic goal is to get know with robust watermarking techniques, especially for medical image marking.

If we look at the types of algorithms, some of the most commonly used (and robust) are Cox, Corvi and Xia algorithm, but all of them rely on watermarks generated as Gauss series of N pseudo-random numbers.

My question is, what if we want to embed custom watermark, like patient information (patient name, date of scan/record, etc.), or to use custom watermark as proof of ownership? Which (robust) algorithms and approaches are used in that case?

Thanks in advance.


texts: Use of gauspuls command in matlab, unable to broad signal at frequency spectrum

texts:
I have a doubt in the gauspuls command in matlab, i entered a fbr(fractional bandwidth of 1.5,ie. bandwidth=  0.75Mhz ) but as i plot the signal in the frequency spectrum using fft. i dont see the signal being broad enough. Why is that. the code is shown below: the below code also has general parameters as follows: dt comes out to be 0.0088us.

%% General parameters %%%

dx=0.05; % mm
Vmax=4.0; % mm/us
alpha=0.99;
dt=alpha*dx/(sqrt(2)*Vmax);% µs
f0=0.5; % central frequency, MHz
t0=1.5; % pulse center time
bndwdth=1.5; % pulse -6dB bandwidth
duration=2*t0; % signal length

timebase=(0:round(duration/dt)-1)'*dt;
[signalI,signalQ,signalE]=gauspuls(timebase-t0,f0,bndwdth);
signal=signalQ;
signal=signal/max(signal);

figure(2)
plot(timebase,signal,'.-')
title('source signal')
xlabel('time (µs)')

fs=1e6;
figure(3)
n=length(signal);
f = fs*(0:n)/n;
Y=fft(signal);
Y=Y';
plot(f,abs(Y))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

Shown below are the plots:


texts: Career Question: As a person who likes math is signal processing a career that's worth it?

texts:
I like high level math and I love deep theory based math.
In industry will I use high level math and theory to implement my own algorithms?
I love high level math but I do not love probability (applications) that much.
Will specializing in signal processing be a red flag?

How does the work of a day to day signal processing engineer in industry is like?


texts: How should you decide on the number of sweeps to be used in FMCW radar?

texts: I am designing an FMCW radar in MATLAB and I would like to know on how should I decide on the number of linear sweeps in FMCW .Does the number of sweeps affect the performance?


texts:
Say I have a PAM signal as $x(t) = \sum_n a_n h(t-nT)$ where $\{a_n\} \in \{\pm 1\}$ are equiprobable random binary bits and $h(t)$ is bandlimited to $[-1/(2T),1/(2T)]$. ($x(t)$ is random process and cyclostationary).

My first question is that can I compute the energy of $x(t)$ as:
$$E_x = E\left\{(x(t))^2\right\} = \int_{-\infty}^\infty \lvert X(f)\rvert^2 df$$
If yes, then my next question is where I get stuck in the solution:
\begin{align}
X(f) &= \int_{-\infty}^{\infty} \sum_n a_n h(t-nT) e^{-j2 \pi ft} dt\\
&= \sum_n a_n \int_{-\infty}^{\infty} h(t-nT)e^{-j2 \pi ft} dt\\
&= \sum_n a_n  H(f) e^{-j2 \pi fnT}\\
&= H(f)A(e^{-j2 \pi fT})
\end{align}
Therefore,
$$E_x = \int_{-\infty}^\infty \lvert H(f)\rvert^2\lvert A(e^{-j2 \pi fT})\rvert^2 df$$

I know the answer is $E_x = \frac{1}{T} E_h E_a$ but I cannot make the connection.


texts: Why we add imaginary part in inverse fourier transform since the time domain signal has only real value?

texts:
The equation of inverse Fourier transform is the following:
$$ f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} d\omega $$
$$ f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) \big[\cos(\omega t) + i\sin(\omega t)\big] d\omega $$
Now, if I have only real value in the time domain $f(t)$, then the imaginary part of $\int e^{i\omega t}$ must go to zero. In our signal processing course, I always had real valued time domain function. Then why am I adding this imaginary part of $\int e^{i\omega t}$ in inverse Fourier transform since it always goes to zero?


texts: How can a signal have two maximum frequency components?

texts: $x(t)$ can be exactly reconstructed from its samples at $\omega_s = 10 \textrm{ rad/sec}$. My conclusion is that the maximum frequency component in $x(t)$ is $5\textrm{ rad/sec}$. But I'm being told that in fact $x(t)$ has maximum frequency component $7 \textrm{ rad/sec}$. How can this be possible?


texts: Kalman Filter | Why Shift the Belief Distribution

texts: I am looking at the video lecture on Kalman Filter. I got some understanding that if the robot moves that we also shift the enviornment belief model accordingly (the same mean amount the robot moved). I am quite confused that why do we move the belief model of the enviornment. The enviornment itself is stationary so if the robot moves the door will remain in place so their distribution should not be shifted. In the following snapshot, the doors are our measurements and we form the belief of where the robot is with respect to the doors. When the robot moves the doors stay stationary. Why do we move the posterior model of as shown in the image below?


texts:
Is there any current ideas about IIR filter algorithm parallelization with thread count greater than filter order?

Ideally, parallelization method should load GPU efficiently.


texts: How to determine load cell response to arbitrary temperature changes with known step response(s)

texts:
We have a loadcell, the output of which is (naturally) temperature dependent.
Under static load conditions, we have observed a dynamic response to temperature changes. Therefore, we did a step response experiment by placing the device in a -20 degree (C) freezer and then transfering it to a 40 degree (C) climate chamber.

What we see is that the temperature of the loadcell has a first order response and the loadcell output has a second degree response (under damped).

What we want to do now is to correct the temperature dynamics induced loadcell output error by using the temperature signal.

How can we use the results we obtained from the step response to correct the weight output when the temperature changes in real time, given temperature and loadcell output signals?


texts: Non-zero DFT components where zero is expected?

texts:
I am implementing DFT in Octave. Here's my code:
function [real_comp, imag_comp] = mydft (samples)
    N = columns(samples);
    for m = 1:N
        real_accum = 0;
        imag_accum = 0;
        for n = 1:N
            real_factor = cos((2 * pi * (n-1) * (m-1)) / N);
            # imag_factor = -I * sin((2 * pi * (n-1) * (m-1)) / N);
            imag_factor = sin((2 * pi * (n-1) * (m-1)) / N);
            real_accum += real_factor * samples(n);
            imag_accum += imag_factor * samples(n);
        endfor
        real_comp(m) = real_acuum;
        imag_comp(m) = imag_accum;
    endfor
endfunction

My goal is to compute the magnitude at each DFT bin and plot it. I am passing in a list of real-only samples and storing the real and imaginary components in different variables. Here's the sample vector that I am using:
samples = [0.3535, 0.3535, 0.6464, 1.0607, 0.3535, -1.0607, -1.3535, -0.3535]

The samples are carefully chosen to prevent spectral leakage at fs = 8KHz. The graph of the DFT magnitude is exactly what I expect. There's nothing wrong with that. But what I don't get is that some components in my 8 point DFT which are supposed be zero, are not exactly zero but rather a very small, close-to-zero value. The built-in fft() function, on the other hand, has zeroed out components where I expect.
For example, I print the elements of output vector produced by my DFT and octave's FFt using the code below:
fft_out = fft(samples);
[real, imag] = mydft(samples);
dft_out = real + (-I * imag);
for i=1:columns(samples)
    disp(fft_out(i));
    disp(dft_out(i));
    printf("\n");
endfor

The point-to-point comparison is given in the table below:

bin_index
fft_out
dft_out

1
-1.0000e-04
-1.0000e-04

2
0.00000 - 3.99988i
4.7184e-16 - 3.9999e+00i

3
1.4141 + 1.4144i
1.4141 + 1.4144i

4
0.0000e+00 - 8.0820e-05i
1.3878e-16 - 8.0820e-05i

5
-1.0000e-04
-1.0000e-04 - 1.4350e-16i

6
0.0000e+00 + 8.0820e-05i
2.3037e-15 + 8.0820e-05i

7
1.4141 - 1.4144i
1.4141 - 1.4144i

8
0.00000 + 3.99988i
-2.0817e-15 + 3.9999e+00i

Now, I wouldn't care too much about the tiny difference but if I compute the phase angles of fft_out and dft_out,
disp(arg(fft_out))
disp(arg(dft_out))

I get the sign flipped in some cases:
fft_out: 3.14159  -1.57080   0.78550  -1.57080   3.14159   1.57080  -0.78550   1.57080
dft_out:-3.14159  -1.57080   0.78550  -1.57080  -3.14159   1.57080  -0.78550   1.57080

So, is there something wrong with my DFT implementation? Is it just a rounding or precision issue?


texts: Phase response in a filter for an Audio Application

texts: Should I be concerned with the group phase delay which a filter produces? Especially in an audio application? Should I aim at using a linear phase delay filter for example?


texts: What are the differences between the DWT and the MODWT?

texts:
What are the differences between the DWT (Discrete Wavelet Transform), which is the most classical algorithm and the Maximum Overlap Discrete Wavelet Trasnform (MODWT)?

Both these algorithm are available in MATLAB. As far as I know MODTW is more suite for multi-resolution analysis, whereas DWT is more useful for signal de-noising and compression.

Thanks for your clarification.


texts: how to extract telephone audio only from an mp3

texts:
I have mp3s of radio talk shows in which the host interviews people on the phone. The host uses a microphone, and the interviewee uses either a landline phone or a cell phone.

I want to save just the telephone audio to a different file. Is that feasible?


texts: Get Spectrogram from Filtered Power Spectrum

texts:
My objective is to get the higher resolution of spectrogram on the high-frequency area (2000 Hz - 5000 Hz) for a section of speech audio. I know that we typically apply Mel-scaled filter bank to get better resolution on lower frequency for speech recognition related tasks. However, in my current work, I want to construct a set of inverse Mel-scaled filter banks to achieve my objective. I have successfully constructed the inverse Mel-filter and applied it to the power spectrum to get inverse MFCC. So I'm a bit lost on this point. How should I get the spectrogram after I applied the inverse Mel-scaled filter to the power spectrum?
Much Thanks,
John


texts: How to Detect a Inhomogeneity Region in Image

texts:
I have an image that includes object and background. However, the object appears some inhomogeneity region due to illumination. My work is that how to detect inhomogeneity region. Which is feature can represent it? Example: low inhomogeneity-high/low feature value....Thank you so much

or


texts: Result of inverse FFT is sometimes shifted in real space

texts:
I am using the Numpy fft2, ifft2, and related functions and I am sometimes running into a strange situation where the output after performing the inverse Fourier transform the results are shifted. I know about fftshift and ifftshift, however my understanding of those is that they are used after fft2 and before ifft2 if it is needed for the DC component to be in the middle of the result instead of at one end. In particular, they are not supposed to be needed in real-space.

One example is in a demonstration of the convolution theorem. The following two things should be identical:

$$F^{-1}(F(im) \cdot circle)$$
$$im \otimes F^{-1}(circle)$$

where $F(x)$ is the Fourier transform, $F^{-1}(x)$ is the inverse transform, $\cdot$ is pointwise multiplication and $\otimes$ is convolution, $im$ is the image, and $circle$ is a simple, ideal, low-pass filter. The code I use to compare these is as follows:

# The fourier duck is a 256x256 8-bit grayscale image
im = skimage.io.imread('fourier_duck.png') / 255

# Create the ideal low-pass filter
x,y = numpy.meshgrid(range(-128,128), range(-128,128))
circle  = numpy.fft.ifftshift(x*x + y*y < 15*15)

# Perform filtering in Fourier space
f_im = numpy.fft.fft2(im)
f_im_circle = f_im * circle
finv_f_im_circle = numpy.real_if_close(numpy.fft.ifft2(f_im_circle), 10000)

# Perform filtering in real space
finv_circle = numpy.real_if_close(numpy.fft.ifft2(circle))
conv_im_finv_circle = scipy.ndimage.convolve(im, finv_circle, mode='wrap') # WARNING: takes lots of memory

The two outputs, finv_f_im_circle and conv_im_finv_circle should be quite similar but are quite different:

If the final, real-space, version is put through numpy.fft.fftshift it looks how I would expect (even though the plots don't show it, the value ranges are essentially the same as well, within FP error). The problem seems to crop up with the inverse Fourier transform of the circle which produces a kernel that also looks shifted.

I have tried doing shifts in other places and looking at the imaginary parts all to no available.

Besides the simple example I show above this has also occured with other filters designed in Fourier space, sometimes when doing inverse filtering, and even with filters designed in real space but applied in Fourier space. I can't figure out the pattern of when to shift the final result or when not to.

Ultimately my question is two fold:

Why is this happening?
How do I fix it? (which may just be how do I predict when I need to use the extra real-space numpy.fft.fftshift)


texts:
I have two related questions on this topic.  First, about how wide should the filter be around the carrier frequency?  What is used to determine this?

Second, in Python, I must use 'numtaps' when implementing a filter.  What exactly are numtaps?  How can I determine what to set this value to?  Setting it to about the carrier frequency seems to make decoding more successful.


texts: Estimating Convolution Input Under the Assumption of Sparsity and Constant Non Zero Values Using Compressive Sensing Approach

texts:
I was wondering about if there is compressive sensing algorithm to estimate the sparse vector where the number of non-zeros values and amplitude of every non-zeros value are known. For example, assume we have a vector $x$ whose length is $N$x$1$ with only $N/2$ equal non-zeros known values but unknown location of those values. That vector  is convolved with a channel $h$ resulting a vector $y$. it means that:

$y = h⊗x$, where $⊗$ is the convolution operation.

Is it possible to use compressive sensing to estimate the locations of non-zeros values in $x$ based on $y$.

Here is the code of an example where the length of vector $x$ is $32$ and channel $h$ = 16:

clear all; clc;

%%%% Build the sparse vector
X = hadamard(32);
X2 = randi([1 length(X)-1], 1);

x = X(1,:) + X(X2+1,:);         % Here the built sparse vector
x = x / max(x);                 % To make the sparse vector either one or zeros

h = randn(1,16);                % channel
y = conv(x,h);
y = y(1:end-length(h)+1);       % To remove the delay of convolution

Thank you