texts:
I am working on a video motion feature implementation. I am implementing a work in a research paper. What I want to implement is this:

  The frames contained in a video clip can be thought of as a volume obtained by considering all the frames in time. This volume can be decomposed into a set of two 2-D temporal slices I(x,t)and I(y,t).

What's the methodology of implementing this i.e slicing the volume of frames in time?
Relevant link of paper: http://www.cs.cmu.edu/~yaser/RasheedSheikhShah_CVST_2005.pdf


texts: Time delay analysis in different binary sequences

texts:
I am extracting a property of a signal which I represent with either a 0 for absent or 1 for present. This property is extracted for each data point in a time series and in the end I am left with a binary sequence for each signal. These sequences are somewhat similar, since their origin is related, but they are different signals (not just a delayed version of each other).
I am currently using cross-correlation to estimate a "delay" between them, but I just read that cross-correlation is mostly suited for delayed versions of the same signal. I would like to know what alternatives I could use.

You can get an example of the data here:

https://www.dropbox.com/s/64uci3juok0qrd2/states.mat


texts: Detecting a three-tone signal in real-time

texts:
I'm new to digital signal processing, and I've been looking into FFT, filtering, Matlab etc. What I'd like is to build a program which will detect when a certain BEEP is transmitted through a trunked voice channel (it signifies an impending critical message) being received via SDR. Through FFT I have ascertained it consists of three frequencies (500Hz, 600Hz and 700Hz). In essence -- somewhat like a phone call with DTMF -- I'd like to monitor audio for the presence of this beeeeep tone and have the software recognize when this is occurring.

My concern is that voice transmissions will cause these three frequencies to peak, making it difficult to determine whether it is due to the alert beep or somebody speaking.

I believe this might be accomplished through some means of measuring the length of the beeps, however this may be a misunderstanding through my inexperience in the field.

I'd appreciate any advice, and if my question isn't directly answerable, please push me in the right direction.

Edit: Here is a temporal frequency analysis of the series from WavePad


texts: Trying to convert Raw image(Bayer format) to jpg/png image in Python

texts:
I am using the Sony Low-light image dataset from the paper Learning-to-see-in-the-Dark. I am trying to convert Raw image (bayer pattern in arw format) into a jpg/png. I used an online tool that gives a very dark jpg for a short exposure image and I have also used Python rawpy library to read and process the raw image. The resulting image from rawpy is shown below:

Rawpy gives a brightened image for a low-light image, which is not desired. I was expecting a dark rgb image. How can I get a dark rgb image in jpg format using rawpy?


texts: Lagrange Multipliers Optimization - Complex Functions

texts:
I want to optimize a sum of complex frequency responses (complex sum is SPL, sound pressure level as a complex function of distance) of a number of loudspeakers, say 10. There is a reference, desired SPL curve given (objective function), which represents the magnitude of the sum of 10 frequency responses. I have read a number of papers how to do optimization using Lagrange Multipliers, however, did not find many papers (only one) about optimization when the functions are complex (Re and Im parts).

Any help or pointing to a good references, examples would be appreciated.


texts: Savitzky-Golay Filter Coefficients and Wikipedia

texts:
I am trying to implement a Savitzky-Golay filter by following this Wikipedia page and in the first instance I have written the following Octave software code to create the convolution coefficients

clear all

m = input( 'Enter number of bars (must be an odd number): ' ) ;

% check no_bars is odd
if mod( m , 2 ) == 0 % is even, report error

fprintf( '\nYou have entered an even number: try again.\n\n' ) ;

else % is odd

% create z
z = ( ( 1 - m ) / 2 : 1 : ( m - 1 ) / 2 )' ;

% create Jacobian matrix J for a cubic polynomial
J = ones( m , 4 ) ;
J( : , 2 ) = z ;
J( : , 3 ) = z .* z ;
J( : , 4 ) = z .* z .* z ;

% create convolution matrix C & show in terminal
C = inv( J' * J ) * J'

end

which gives a matrix output where the first line are the coefficients for the smooth, second line are coefficients for the 1st derivative, third line are coefficients for the 2nd derivative etc. However, when I check this program's output with the selected coefficients given in the Wikipedia page appendix there is disagreement.

To do this check I take the normalisation figures given in the tables, e.g. 35 for smoothing coefficients for a window of length 5 for a cubic fit, and multiply the matrix by this number, which will then transform the relevant row of the matrix into the coefficients given in the tables. Following this procedure I see that for the actual smoothing coefficients and the coefficients for the 1st derivative for window lengths 5, 7 and 9 for a cubic fit there is agreement. However, for the 2nd and 3rd derivatives there is no agreement. In fact to get agreement the normalisation figures in the table would need to be changed thus:

window length 5 - 2nd deriv = 14, 3rd deriv = 12
window length 7 - 2nd deriv = 84, 3rd deriv = 36
window length 9 - 2nd deriv = 924, 3rd deriv = 1188

Therefore I am now uncertain whether I have understood the Wikipedia page and have correctly implemented the above code, or whether in fact the Wikipedia page is wrong in this regard. Can anyone point out where the problem is?


texts: How to create a distortion compensation filter

texts:
I am creating a channelizer that consists of a complex mixer, CIC decimator, and a FIR compensation/decimation filter.  The final FIR filter may, if it matters, be implemented as more than one filter.

My question is, how do I design a filter such that it compensates for the very non-flat frequency response of the CIC filter?  Do you create the frequency response that you want by calculating the reciprocal of the CIC's response and then run it through an inverse FFT to get the impulse response?

As you can see, though my particular problem pertains to CIC filters, my question is really about how you create any kind of distortion compensation filter.

Thanks for your time.


texts: HHT and how to plot the Hilbert spectrum

texts:
NOTE:
The following is based on an unanswered question posted here:
How to plot the Hilbert Spectrum in Hilbert-Huang transform?
, based on the usage of the "plot_hht" from here:
(https://www.mathworks.com/matlabcentral/fileexchange/19681-hilbert-huang-transform)

From what I understand the "plot_hht" does,
after applying EMD to a signal, the IMF's of the signal remain.
The Hilbert Transform is then applied to each IMF, and the resulting phase angles are somehow used to help find the instantaneous frequency (with respect to time).

Is it possible to represent these instantaneous frequencies vs time, with a third dimension, its strength/amplitude?

Thanks in advance.

EDIT P.S.: For clarity, the desired result is a graph that looks like this:
,
as taken from an article titled "Detecting position dependent tremor with the Empirical mode decomposition".


texts: Why my result from FFT which is negative frequencies are different from positive frequencies

texts:
I am trying to analyze a signal from monkey brain by applying FFT on the signal.
By using numpy.fft.fft function to perform Fourier Transform on raw signal, I got the power of frequencies in the bellow figure.

My question is why the power of negative frequencies are different from positive frequencies.

Please help me to give your opinions if you have ideas about it.

Thank for helps.


texts: PortAudio: sound fade in/out by scaling sample data

texts:
I've done already some fade in or fade out effect in portaudio, but that is not working properly.

My callback function:

int SoundEngine::patestCallback(const void *inputBuffer, void *outputBuffer, unsigned long framesPerBuffer, const PaStreamCallbackTimeInfo *timeInfo, PaStreamCallbackFlags statusFlags, void *userData)
{
    paData *callData = (paData*)userData;
    float *out = (float*)outputBuffer;
    float sample;
    unsigned long i;

    (void) timeInfo; /* Prevent unused variable warnings. */
    (void) statusFlags;
    (void) inputBuffer;

    for( i=0; i < framesPerBuffer; i++ )
    {
        sample = callData->sine[callData->phase++];
        *out++ = sample;

        if( callData->phase >= SE_TABLE_SIZE ) callData->phase -= SE_TABLE_SIZE;
    }
    return paContinue;
}

I've discovered that for whole stream it's easy to make fade out effect, by applying something like this:

if(callData->fadeFlag == true)
        {
            callData->toneCountdown--;
            sample *= (-cos(((float)(callData->toneCountdown - (callData->toneLength - callData->fadeOut)) / (float)callData->fadeOut) * M_PI) + 1.) * .5;
        }

But this is not working for me. If I have more than one sound at the moment, there are fading all sounds. At the moment I'm using callback function as I wrote above. callData->sine at the beginning contains only 0 to not produce any sound. Sound I'm adding by calling a member function playSound(). I am doing this in single thread - so I can use a lot of these functions at the same time - making a multiple sounds. How is working function playSound? There it is:

void SoundEngine::playSound(int instrument, int freq, int duration)
{
    if(instrument == 0)
    {
        for(int j=0; j<20; j++)
        {
            for(int i=0; i<SE_TABLE_SIZE; i++)
            {
                data.sine[i] += 0.01*sin(2 * M_PI * data.soundPitch[freq] * ((double)i/(double)SAMPLE_RATE));
            }
        }
        Pa_Sleep(duration);
        for(int j=0; j<20; j++)
        {
            for(int i=0; i<SE_TABLE_SIZE; i++)
            {
                data.sine[i] -= 0.01*sin(2 * M_PI * data.soundPitch[freq] * ((double)i/(double)SAMPLE_RATE));
            }
        }
    }
}

You are probably asking now: "What the hell that man done here?". I know that it's not a masterpiece of coding, but I'm newbie in sound programming, so I'm trying my best to learn it. For now I want just make it work as it should work.

Function playSound() is working great, but fade in/out is making some noise during fadeing. Fade effect is made by for(int j=0; j<20; j++) that's gradually volume up my sound. I think that modyfing data.sinedata couse those noises. When I slow down fade effect - it works perfectly, but too slow. I need fast fade in/out effect to eliminate noises at start and end of sound.

Can someone help me to apply fade in/out effect to my example?


texts: Peak location estimation for a triangle wave in AGWN?

texts:
Is there a more optimal peak location estimator specifically for triangle waves in AGWN, compared with peak estimation for generic waveforms?

For a "smooth" waveform in noise, a common peak locator method might involve a linear-phase low-pass filter or an overdetermined polynomial curve fit (by regression, et.al.), followed by a maxima/minima search.  But if the signal is known to be a sharp triangle wave (which may imply that any sampling included/aliased some non-filtered frequency content above Fs/2, if that is relevant), or portion thereof, and not a low order polynomial or other smooth waveform, can this information be used to improve the triangle peak location estimation?

If so, how?

Added: Less than 1 period of the triangle wave may be available.  But assume at least one extrema is in the data window.


texts: Alternative to BIBO stability of a system

texts:
In DSP textbooks a system is stable in the BIBO (Bounded-Input, Bounded Output) sense if and only if every bounded input sequence produces a bounded output sequence. After stating this definition stable always means BIBO-stable.

I'd like to know if there exists other form of stability for a system?


texts: What is bit rate? How is it related to bandwidth?

texts: I am very confused about when to use bandwidth as bit rate and where not. Are these terms one and the same?


texts: How can I call the area under an impulse/spike?

texts: I am looking for an appropriate term to express the area under a spike or impulse for time series signals. The term "Briefness" is probably not appropriate since the amplitude should be taken into account


texts: Convolution engine for FIR and IIR equal?

texts:
I have a beginners question. I just installed Audiovero Acourate. This software create FIR filters for room correction of my loud speakers. The filters are .wav files which then should be used by a convolution engine, like jconvolver or brutefir on Linux.

Now, I decided to try to setup a convolution engine on Windows instead, since there are several pain points (software and driver / hardware) on Linux. I now use Equalizer APO on Windows for the convolution.

My question is: Is there any difference between convolution engines for IIR and FIR filters? And is there a difference in sound quality using either convolution engine?

I could not find much info about the convolution engine within Equalizer APO. I‘m just interested in best audio quality.


texts: Extracting top & bottom fields from a YUV 4:2:0 Frame

texts:
I am wondering what would be the right approach to extract top and bottom fields from a YUV 4:2:0 frame.

In the YUV 4:2:0 representation, Luma and Chroma samples are located as shown below:

Y Y Y Y Y Y Y Y
c     c     c    c
Y Y Y Y Y Y Y Y

Y Y Y Y Y Y Y Y
c     c     c    c
Y Y Y Y Y Y Y Y

Hence, one needs to be careful while extracting chroma top and bottom fields (due to the chroma offset problem). I have tried a few methods (using asymmetric filters and such), but none have been satisfactory.

Are there some good resources (papers) on this problem?


texts: Would subtracting the phases of two images be a superior difference metric than subtracting the images directly?

texts:
I'm hoping someone can sanity check this idea as I am admittedly a bit of a noob when it comes to working with FFT's.

Say I have two (registered) images of two of the same object and want to use one as a baseline to check for differences in quality control.  My initial naive approach is to just subtract the two images directly and then whatever remains can be treated as defects.  However, this method is subject to error in the presence of illumination differences.

I'm thinking that I would be better served if I take the FFT of the images, and subtract only the phase information, and then use the IFFT of that result as the defect map, and that this should hopefully eliminate false positives due to lighting.

Does this seem like a reasonable assumption or is there some kind of detail that I'm overlooking?  Thanks for any advice you can offer!


texts: Second moment properties (covariance) of lag window estimator

texts:
We know that the lag window estimator of the spectral density is defined as:
$$S^{(lw)}(f)=\int_{-f_{(N)}}^{f_{(N)}}W_m(f-\phi)\hat{S}^{p}(\phi)d\phi=\sum\limits_{\tau=-(N-1)}^{\tau=N-1}w_{\tau,m}s^{(p)}_{\tau}e^{-2\pi f\tau}$$
Where $W_m(f)$ (the spectral window) is the Fourier transform of $w_{\tau,m}$( the lag window), and $\hat{S}^{p}$ (The periodogram ) is the Fourier transform of $s^{(p)}_{\tau}$ (the covariance function).
 The covariance between  $S^{(lw)}(f)$ and $S^{(lw)}(f')$ is approximatly $0$ if $|f-f'|$ is larger than the width (or bandwidth) of the lag window $W_m$.  But if we write the following:
$$cov\{S^{(lw)}(f),S^{(lw)}(f')\}=E\{S^{(lw)}(f)S^{(lw)}(f')\}- S(f)S(f')$$
$S(f)$ being the true spectral density at frequency $f$. We can write then:
$$E\{S^{(lw)}(f)S^{(lw)}(f')\} =\frac{1}{N}\sum\limits_ {f,f'}W_m(f)W_m(f')E\{\hat{S}^{(p)}(f)\hat{S}^{(p)}(f')\}$$
$$=\frac{1}{N}\sum\limits_ {f,f'}W_m(f)W_m(f')E\{\hat{S}^{(p)}(f)\}E\{\hat{S}^{(p)}(f')\}=S(f)S(f')$$
because we know that the periodogram ordinates are uncorrelated for $|f-f'|> \frac{1}{N}$. From what I have just written the ordinates of the lag window estimator are uncorrelated onlu if the peridodogram ordinates are uncorrelated ($|f-f'|> \frac{1}{N}$), which is inconsistent with the result that we already know. Where am I mistaken?


texts: What is the purpose of the log when computing the MFCC?

texts:
The steps of computing the Mel-Frequency Cepstrum Coefficients (MFCC) are:

Frame blocking -> Windowing-> abs(DFT) -> Mel filter bank-> Sum coefficients for each filter-> Logarithm -> DCT

But what is the purpose of the logarithm step?


texts: Signal Timing Recovery, Loop Bandwidth and SNR

texts:
As explained in the answer here, for the symbol timing recovery loop we need to keep loop bandwidth sufficiently below the modulation bandwidth of the signal so as to avoid the filtering of the signal itself.
But when we send a random sequence of bits the corresponding baseband signal will contain frequencies from $-\frac{f_s}{2}$ to $+\frac{f_s}{2}$, where $f_s$ is the symbol rate.
Since the signal has spectrum even at low frequencies, how can we ensure that the signal will not be filtered by the timing recovery loop?
I am asking because of the following reason:
1. I generate the PSK modulated data and add phase noise on it

dataMod = pskMod(dataVec, 8, 'gray');
tx_out = dataMod.*exp(1i*ph_out);

Here, the ph_out is the phase jitter signal and dataMod is the PSK modulated signal of the dataVec symbol stream.
2. I pass the phase of the data to the carrier recovery loop:

pll_out = pll(angle(tx_out), 0.005, 1e-6, 1e-3, 1);
pll_out1 = pll(ph_out, 0.005, 1e-6, 4e-3, 1);
pll_out2 = pll(angle(dataMod), 0.005, 1e-6, 1e-3, 1);

I create three different outputs, pll_out is the loop output with the complete signal, pll_out1 is the loop output with only the phase noise as input and pll_out2 is the output with only the signal (without noise) as the input.
I see the following at the output of the tracking loop:
1. The power spectral densities for input and output phase noise:

From this, I can see that the low frequency phase noise is getting filtered as expected.

The power spectral densities of input and output signal:

From this, it is seen that the signal power at low frequencies is filtered. But because of this filtering additional "noise" or "distortion" (whatever is the correct term) causes the following PSD at the output (yellow):

The frequency roll-off is explained very well by Dan here. But somehow now, the roll-off is not the problem for me.
3. The total noise at the PLL loop output is thus:

So, although the signal distortion/noise does not affect me at high frequencies (and its roll-off doesn't matter anymore) but my low frequency noise gets degraded.
So, how can we ensure that most of the signal power is outside the PLL bandwidth?


texts: Good way to detect pulse with known width with background noise?

texts:
I'm dealing with a light sensor tracking a bursting LED. Unfortunately signal is weak and there are other light sources in environment. I've used optical filters on sensor but there is still some noise.

The LED burst with a fixed duration but the distance and time vary. So I'm interested in finding those 5 pixels wide pulses (the height of pulse may be different). The noise is often a square wave with a much lower frequency.

What algorithm can be used to detect those pulse?

Thanks.


texts: Understanding a Bainter Circuit (chebychev filter)

texts:
As part of a class in Signal Processing, I am building a 3rd order Chebychev band reject filter. We implement this by using three cascading Bainter circuits. Although not part of the class, I have a question about the gain of the Bainter circuit.

I am trying to write a script which will automate the component selection using the corner frequencies and maximum overall gain as design rules, but have some issue with calculating the overall gain.

To calculate the overall gain of a Bainter stage, would I simply work out the individual gains of the three op-amp sections? The overall gain would then be the product of the three individual gains?


texts: What information does the abs(normalized_image - normalized_quantised_image) carry?

texts:
Given a grayscale image X2

And its (under)quantised version with 3 bits depth X3

What information does the abs(X2_normalized - X3_normalized) carry?As an example the following image X5 generated with the code below

X2norm= X2-min(X2(:)); X2norm = X2norm / max(X2norm(:));
X3norm= X3-min(X3(:)); X3norm = X3norm / max(X3norm(:));
X5 = abs(X2norm - X3norm);
image5=figure;imagesc(X5);title('image 5');
saveas(image5,'image5.jpg');


texts:
I am trying to apply Haralick textures to a SAR image (float32). As far as I know, the image first needs to be quantized to a reasonable bit depth prior to calculate the co-ocurrence matrix. In the original paper, Haralick proposes using a "equal probability quantizer". As far as I know is the algorithm employed by Mathlab "imquantize". I have been looking for a pre-existing python/numpy implementation (no success).
I also have reading pre-existing java code from ESA SNAP. As far I know histogram calculations would be straightforward with numpy, but if I understand correctly, every pixel is assigned  a "level/histogram bin" using bisection search, so I do not see any straightforward way to properly vectorize the code.
EDITED after this point ===================
I have made simple implementation in Python + Numpy based on the pre-existing java code from ESA SNAP, and the kind advice of Marcus Müller.
import numpy as np

def eqProbQuant(image, levels=32):
    # Sort the pixels by value
    sorted_image = np.sort(image)
    # Get the pixel count
    pixel_count = sorted_image.shape[0]
    # Get the number of pixels per bin
    samples_per_bin = int(pixel_count/levels)
    # Get locations where the bin would change
    edge_samples = np.arange(levels+1) * samples_per_bin
    # Get the values at those locations (bin edges)
    bin_edges = sorted_image[edge_samples]
    # Use the values to apply quantization
    quantized = np.digitize(image, bin_edges)

    return quantized

If possible, I would like to know if this code would apply the algorithm correctly (I have very little experience with signal processing)


texts: Confused about the Fourier Transform of $e^{at}u(t)$

texts:
This is the problem at hand:

I'm unaware of why we didn't have to say anything about $\omega$ like that it should be also greater than $0$, I know it's variable...but it's multiplied by $t$ whatsoever.
Also, if $a$ is complex, I don't understand why do we require that $Re \{a\} >0$ only, without mentioning anything about the imaginary part (it's mentioned further down in the book).
Finally, why isn't it $Re\{a\}>=0$? is it because that may reduce to the Heaviside step function if the imaginary part of $a$ is also 0?
Thank you.


texts: Brief introduction to wireless communication and communication terminologies

texts:
In the Question Design of equalizer for wireless communication the diagram of a basic communication system is presented with informative points and discussions. My Questions are:

Considering the channel modeled as a FIR filter having a single input and single output. If the impulse response (channel coefficients) take values in the complex domain then am I looking at a Rayleigh Fading channel?
What constitutes a multi-path channel? Would multi-path channel be typically not a single input and single output?
In the diagram, the input to the channel is Tx. Is Tx digital? I mean can it take bits or multiple level bits? Suppose, we want to send an image or a real-valued time series through a communication system, then is it transmitted in real values or digital in wireless communication?
At which stage is digital modulation performed? Then the output of the modulation should be digital. So, is this digital output the Tx that is sent to a channel?
What is the difference between modulation and encoding? When we are doing encoding, we are basically transforming the information from one representation to another. Isn't this what modulation also does?


texts: Autocorrelation and FFT : avoid zero-padding

texts:
Some posts evocate the computation of convolution or cross-correlation using FFT and zero-padding the temporal signal.

I want to compute the autocorrelation of a 3D array using FFTW (2 dimensions are periodic). The file is very large (dimension 256*256*49995, double precision, 25GB).

1 : Is padding necessary for periodic dimensions?

2 : Is it possible to avoid zero-padding to limit memory usage? (fftw++ seems to answer yes)

3 : From the 3d FFT of the non-padded signal, is it possible to compute the cross-correlation using a Discrete Cosine Transform?

Cross-correlation = DCT[ FFT(signal)*Conjugate(FFT(signal)) ]

EDIT - 4 : Is it mandatory to use a mask?

Thanks


texts: How to determine cut-off frequency using power spectral density?

texts:
I am new to DSP and my question is a very simple one. It is about determining the right cut-off frequency. I calculate the power spectral density using FFT by Welch's algorithm. However I am confused as to what should be the right cut-off frequency.

UPDATE

Can anyone tell me what is the cut-off frequency in this case:


texts: Does Wigner-Ville distribution provide the highest resolution?

texts:
For linear time-frequency representation, the resolution is limited by Gabor limit. However, if interference is allowed, the resolution may be higher. What is the possible resolution in such a case? In "Electroencephalography Feature Extraction Using High Time-Frequency Resolution Analysis ", they stated

  We  all  know  that  the  Wigner-Ville  distribution  achieves  to  the  lower  bound  of  the
  uncertainty  principle  and  have  the  highest  resolution  in  all  time-frequency  distributions

Why is this statement true?


texts: Non Circulant Translation Using Fourier Transform

texts:
The translation property of Fourier Transform (FT) for a two dimensional image $f$ is as
$$
f(x-x_0, y - y_0) = F(u, v)e^{-j2\pi(ux_0/M+vy_0/N)}
$$

Using this equation, the following code (in Matlab) computes and shows translation of the input image.

im = im2double(imread('cameraman.tif'));
F = fft2(im);
[r, c] = size(F);
x0 = 20;
y0 = 40;

% [u, v] = meshgrid(1:c, 1:r);
[u, v] = meshgrid(-c/2:c/2-1, -r/2:r/2-1);
A = fftshift(exp(-1i * 2 * pi * (u * x0 / c + v * y0 / r )));

Y = F .* A ;
y = ifft2(Y);
imshow(real(y), [])

The output of above code is as follows

but I need output image be as following

My question:
How I can produce The second image using FT?
Thanks.


texts:
For an exercise, we made a noise signal by generating a normally-distributed random vector in MATLAB. The mean is defined to be equal to zero and the variance is set to 0.5.

I know that the power spectral density is just the fourier transform of the variance in frequency domain. In this case, does the PSD just equals to 0.5 since the variance is constant with time?


texts:
Hi I am a little confused on what the notation of the following statement means.

$$ H_{k}(z)= H(W_{4}^{k} z), k = 0,...,3$$

It comes from a question in which I have designed a FIR low-pass filter $H(z)$ and my goal is to implement a DFT filter bank scheme like this:

Exchange $P(z)$ for $H(z)$ and k would correspond to the subscript of P and M in this case is equal to 4

I guess I am confused on how to find $H_{k}(z)$ or what exactly a polyphase filter is


texts:
I was trying to find the period of $\sin(\frac{12\pi n}{5})-\sin(\frac{2\pi n}{5})$, each of the sinusoids has a period of 5 however their difference has a period of 1. It turned out that they're identical in the discrete-time domain.
besides actually solving the equation for reals and concluding integer roots, is there anything that explains this in terms of periodicity?
Thank you.


texts: Is Sum of Absolute Value / $ {L}_{1} $ Norm of Differences Convex?

texts:
I'm not sure how to approach this exercise.
One idea is to derive it w.r.t z, show that there is a min-extremum at $z=f_k$ and then show that for each value from the right and the left of the loss function it is positive which will prove that it is convex.
I never worked with Dirichlet functions, I'm not sure how should I show that it is convex analytically. Can somebody give me a clue/solution path ?


texts:
For practicing spotlight processing, i am looking for I/Q baseband signal from a SAR airborne or spaceborne system. In my research I have only found processed image from :

eoweb geoportal https://eoweb.dlr.de/egp/
Airbus Terrasar X https://www.intelligence-airbusds.com/imagery/constellation/

Could anyone give me a link for raw I/Q signal (and associated metadata)?
Thanks a lot


texts: Do I have to flip my kernel when performing an FFT-based convolution?

texts:
I have a 2D image and I convolve it with a 2D kernel image using FFT. So far I was always using symmetric kernels (e.g., Gaussian with stddev_x = stddev_y). As a result, I never bothered thinking whether I have to flip my kernel image or not because it wouldn't have made any difference. However, now I want to convolve my image using an elliptical Gaussian kernel with stddev_x != stddev_y and an arbitrary angle. Now, do I have to flip my kernel image prior the FFT convolution? Or the flipping is required only when using the usual convolution algorithm and not the FFT-based one?

Thank you.


texts: What's the difference between nonstationary and time varying?

texts:
I am confused by the following statement from the paper Beamforming: A Versatile Approach to Spatial Filtering by Barry D. and Kevin M.:
There are two basic adaptive approach: 1) block adaption ..., 2) continuous adaption.... If a nonstationary environment is anticipated, block adaption can be used, provided that the weights are recomputed periodically....Continuous adaption is usually preferred when statistics are time-varying or (for computational reasons) when the number of adaptive weights M is moderate to large (values of M>50 is not uncommon).
This statement gives the conditions when the two kinds of adaption methods should be used. But the condition NONSTATIONARY and TIME-VARYING, in my opinion, are the same. Can you point out in what situations they are different? Illustrative examples would be greatly appreciated.
Thank you.


texts: Emotional Tone detection using audio on mobile phones

texts:
Which audio features should be used to detect tone from audio? Are there any open source frameworks or implementations which can be ported to mobile phones to process audio?
Alternatively, would it be possible to collect audio features from the phone and then process this on the server so as to detect tone?


texts: How to prove main-lobe side-lobe ratio in hann window is $-32\textrm{ dB}$?

texts:
I found that the value of the main lobe is about $N$, when I try to calculate the value of the sidelobe I got $0$, I assumed that the peak is at  $\theta=\frac{6\pi}{N}$.

How can I find the size of the
side-lobe?


texts: Why use windowing function instead of truncating the signal to integer amount of periods?

texts: According to this source of information if the signal is consisted of integer amount of periods then FFT would work ideally as if the signal was infinite in time. The same source mentions that windowing is not completely eliminating the spectrum leakage but just attenuating it. Therefore, my question is why perform windowing instead of truncating the signal to integer amount of periods, shouldn't the truncation to integer amount of periods give better results?


texts: changing filtering "speed" (cutoff frequency) of an IIR filter knowing only its coefficients

texts:
i got an IIR filter bu I have only the coefficients.
now i'd like to be able to change the filter characteristics of the low pass filter (the "cutoff" frequency) but all i got are a and b coefficients.
i'd like  to be able to set a multiplier that will "scale" the filter by that amount. if i set the multiplier to 0.5 i'd like to filter the same but two times lower or two times slower (50 Hz cutoff instead of 100Hz)
is it possible ? Or is it a lost cause ?
and with Matlab ?
thanks so much
Jeff


texts: regarding autocorrelation and derivation of original signal

texts:
$$r_x(t)=\begin{cases}
2-2t    ,&   0 \le t \le 1\\
2+2t ,& -1 \le t \le 0
\end{cases}$$

So how do I get original $x(t)$ from this?

I tried integrating the equations, and the value I got was $1$.


texts: I have data that is affected by temperature.  How can I filter out the temperature signal considering I have temperature data collected?


texts:
I am going to be straight forward and say this is one of the problem sheets and I am stuck on one of the questions.  I already know the answer as I have been given them but I am trying to work backwards to see how to derive from the answer.  so the question is:

 Find another cosine y(t), whose frequency is as close as possible to that of x(t),
 which when sampled with T=1ms yields the same sample values as x(t).
 Write the equation and sketch y(t) vs. t.

where $x(t) = 10\cos(500\pi t)$

here is the answer

the question is I don't quite understand why $10\cos(1500\pi t)$ would go through those points and not $10\cos(1000\pi t)$.


texts: Using Kalman filter for high-dimensional observation

texts: My observation space is essentially a noisy high-resolution video. In this case, regular Kalman filter doesn't work because the covariance matrix is of size (pixel x pixel) and just doesn't fit in memory. What is the typical solution in this case?


texts: I have data from a power transformer that has a frequency of 60 hz.  Sometimes there are faults in the system and causes a large spike in the reading of the transformer.  What I want to do is run the data into a signal processor to detect when these large spikes occur and either remove the spike value or present the time of occurrence.  how would I be able to detect this spike in Octave?


texts: Notch Filter - Bandwidth vs pole / zero location vs Q factor

texts:
I am reading up about notch filters. I understand the closer the pole is to the zero, the narrower the width gets and thus the bandwidth gets smaller. Is there an equation that links the bandwidth to the location of a pole? I have tried searching for so long and I couldn't get any useful answers.

I want to be able to determine the position of the pole given the required bandwidth and sampling rate and the frequency of the noise. I do not want to do it via trial and error

Thanks


texts: How to differentiate the product signal $f(t)\theta(t)$, where $\theta(t)$ is Heaviside's unit step function?

texts:
What is the derivative (in the engineer's sense) of the causal function $f(t)\theta(t)$, where $\theta$ is the Heaviside unit step function?

I've seen the formula $f'(t)\theta(t)+f(0)\delta(t)$, where $\delta$ is Dirac's delta. This looks like a kind of "product rule": differentiating the product gives $f'\theta+f\theta'$, but $\theta'$ is $\delta$, and $\color{blue}{f(t)\delta(t)=f(0)\delta(t)}$.

If this is right, I don't understand the following argument, from the solutions manual to Oppenheim and Wilsky's Signals and Systems. The solutions manual says the derivative of the function $2e^{-3t}\theta(t-1)$ is

$$-6e^{-3t}\theta(t-1)+\color{red}{2}\delta(t-1)$$

It's the second term I don't understand. Using the "product rule" heuristic, the second term should be $2e^{-3t}\delta(t-1)$, which using the blue formula above gives $\color{red}{2e^{-3}}$ times the delayed delta function, not just twice the delayed delta function.

Is the solutions manual wrong?

(Cross-posted on MSE.)


texts: I went through a paper in random demodulator. It states that If the sampling rate is $R$ ($R$ is less then $W$, $W=$band limit of a signal in Hz), and assume that $W$ is divided by $R$. Then each sample is the sum of $W/R$ consecutive entries of the demodulated signal. Does it mean that the signal is compressed by taking $W/R$ samples at a time?


texts:
I am trying understand LTE channel estimation. I have multiple questions:

How much phase shift does a typical LTE channel have?
Can all that be corrected using pilots?
Where do we typically implement these algorithms? In a processor or FPGA? Because the OFDM system will be in an FPGA or similar platform.


texts: Relation between the matrix trace and the amplitude of each element

texts:
Assume a diagonal matrix $\mathbf X$ whose size $N\times N$  and its diagonal elements are $0.5 + 0.5i$, and the vector $\mathbf p$ of size $N\times 1$ whose elements have similar amplitude.
I have noticed that
$$\operatorname{tr}\bigg\{\left(\mathbf{X^{-1}\cdot p\cdot p^H\cdot \left(X^H\right)^{-1}}\right)\bigg\}\tag{1}$$
is related into the amplitude of $\mathbf p$. where $\operatorname{tr}$ is the trace operator.

For example if $\mathbf p = 1$;
$$\implies\operatorname{tr}\bigg\{\left(\mathbf{X^{-1}\cdot p\cdot p^H\cdot \left(X^H\right)^{-1}}\right)\bigg\} = 2N$$

On the other hand, if $\mathbf p=0.5$;
$$\implies\operatorname{tr}\bigg\{\left(\mathbf{X^{-1}\cdot p\cdot p^H\cdot \left(X^H\right)^{-1}}\right)\bigg\} = \frac N2$$

And so on.

My question, is there a general formula expressing the relationship between the  equation $(1)$  and the amplitude of $\mathbf p$ ?


texts: Does this Lucas-Kanade registration MATLAB script appear to be correct?

texts:
I'm currently writing a MATLAB script to replicate the Lucas-Kanade registration algorithm discussed here: https://cecas.clemson.edu/~stb/klt/lucas_bruce_d_1981_1.pdf
Specifically, I'm attempting to replicate the 2D implementation:

However, I'm struggling to get an output that accurately describes the translation I'm inducing.
First, I define two partial derivative functions "x_derivative" and "y_derivative" which simply take the difference of the input signal $f$.
Second, I use the derivative functions to calculate the X and Y partial derivatives.
Third, I input the matrices into the equation to calculate the X and Y components of the H vector. However, I appear to be getting "bad" results.
For example, I've been inputing a signal $g$ that is a 2 pixel right shift of $f$. However, I'm getting "bad" results. Ideally, I should be getting a vector H = [2, 0].
Any obvious mistakes I'm making?
function [h_x, h_y] = estimate(f, g)

% declare derivative function handles
x_derivative = @(x)[diff(x, 1, 2), zeros(length(x(:, 1)), 1)]; % we append zero column such that output dimension = input dimension
y_derivative = @(x)[diff(x, 1, 1); zeros(1, length(x(1, :)))]; % we append row column such that output dimension = input dimension

f_derivative_x = x_derivative(f);
f_derivative_y = y_derivative(f);

% estimate x component
x_mat_numer = transpose(f_derivative_x) * (g - f);
x_sum_numer = sum(sum(x_mat_numer));

x_mat_denom = transpose(f_derivative_x) * f_derivative_x;
x_sum_denom = sum(sum(x_mat_denom));

h_x = x_sum_numer / x_sum_denom;

% estimate y component
y_mat_numer = transpose(f_derivative_y) * (g - f);
y_sum_numer = sum(sum(y_mat_numer));

y_mat_denom = transpose(f_derivative_y) * f_derivative_y;
y_sum_denom = sum(sum(y_mat_denom));

h_y = y_sum_numer / y_sum_denom;

end


texts:
I'm hoping someone could point me in the direction of some signal processing methods to clean up my data. I'm collecting physiological data from human muscle (lower leg; gastrocnemius) during walking/gait.

I know the sensor is picking up contractions, but the problem is it's so plagued by motion and skin artefacts I can barely do anything with it. Performing an isometric contraction (contracting the muscle without changing its length) shows a good signal (as there are no motion or skin artefacts to contaminate it. However, when collecting real-world data of a subject walking I'm getting huge spikes on the heel strike, amongst other unpleasantries.

The sensor type revolves around a MEMS microphone, so basically it picks up everything. It is un-amplified and only filtered by post processing. The sampling rate is 1kHz and seems to be more than enough. It's a single channel/vector of a continuous signal (i.e. 10 minutes of collection at 1kHz is a 1x600000 vector).

Please help me by suggesting some methods on what I can do. I've been looking at wavelet transforms which look promising, but any help on the matter will be much appreciated. Also look at PCA but having a vector makes it difficult (unless I break my data into windows?). I'm trying to keep it as uncontrolled as possible, so any analysis that I can do on my current collected data would be great!

Can provide some data if it helps in your analysis.

Cheers!


texts: I know that many books and papers talk about the DC offset/DC component of a filter. How do we define the DC offset mathematically, for the case of discrete filters?


texts:
I'm used to filtering using FIR filters with a cutoff and dB attenuations etc.
The issue I have is that I have pulses in my signal which contain similar frequencies as the noise. I'd like to totally remove the noise whilst keeping the pulses intact (no attenuation).
What I'd like to do is have a filter which is applied if the frequency content within the window is below a given magnitude and not if it's above (or really an exponential weighting depending on the magnitude of the frequency present in a given window.)
I can probably make this happen in a fairly hacky way - but is there already a developed area of DSP / mathematics - which deals with this issue? - I rather not reinvent the wheel.
It would be nice if this was fast and I didn't have to FFT each window to check the frequency / magnitude.


texts: How to determine the original size of the FFT used in unknown OFDM symbol

texts:
Given an unknown OFDM symbol consisting of 9144 samples, interested to know the size of the FFT used in the TX/Number of bins.

Down-sampling to 512,256,128,64,32 yields the following.

Do not see how can this information be derived based on this data.
As a result, have tried to analyze the time domain of the downsampled vectors, and detect samples that were removed by the downsampling process.

Seems that the third graph represents the original FFT size – 128, as the sample at 277 has not been witnessed before(fft=64,32), and is witnessed at higher sampling rates.

Wonder if this is correct and whether there is a more accurate approach to determine the FFT size used?


texts: Phase difference calculation for non divisible sample rate to frequency

texts:
I am learning DSP and I am trying to simulate a simple phase difference calculator between two sinusoidal signals using GNU radio.

My flowchart looks like this:

As can be seen on the chart, one of the signals is using a delay block (controlled by the "delay sig" range GUI).

Other interesting points are:

fft_width variable:
Calculated by an expression - int(freq/20)
freq variable:
Controlled by a GUI range with steps of 1.2KHz. Min 1.2KHz, Max 9.6KHz
The skip head and Keep M in N blocks - used together to extract only the FFT bin in interest. To do so the skip head calculates the first N blocks to be skipped via this expression: int((samp_rate/fft_width))+2 . Idea taken from this video: https://youtu.be/GJKbD--MsLM?t=326

The problem:

This flowchart works when the samp_rate and freq are divisible - say samp_rate=10k and freq = 1k. However, for more of a real world scenario I have chosen variables that have lower ratio and are not divisible. In particular: samp_rate=11k and freq=1.2k. As a result of this the phase calculation fluctuates between 180 degree difference - for example 3.42 and 3.42-pi.

The question:

How do I make this flowchart work when the samp_rate and freq are not divisible? Furthermore how do I make it work for samp_rate of 11k, but different frequencies (chosen by the freq range GUI)?

What I tried:

Not much unfortunately, as I cannot find much information about this problem. That said, I read somewhere that clock synchronization via a PLL or Clock Recovery MM might help, but I do not know how to apply the relevant GNU radio blocks to this flowchart. Furthermore, if PLL is used would it not change the phase difference between the signals when trying to synchronize them?

Thank you for your time!

EDIT:
I tried putting a PLL ref output block between the throttle and FFT blocks:

For the parameters I used:

loop bandwidth: pi/200 (documentation states it should be either that or 2*pi/100)

min_freq & max_freq: (2*pi*1.1e3)/samp_rate (documentation states it needs to be in radians per sample)

My goal was to try to make it work for 1.2k frequency, but unfortunately I still face the same problems. How can I fix it?


texts:
I've been learning about digital signal processing for some time, and I understand the principles in the digital world, as well as the conversion of continuous filters in the s-domain to digital filters in the z-domain.

I'd love to learn more about the process of designing analog filters and converting this into circuits that I can try and test. Since I cannot brag about my knowledge in electronics, with the exception of a little experience dabbling with Arduino projects, I don't see how an equation in the s-domain can be translated to a circuit. From what little I have found, most sources go the other way, taking a circuit and analysing it in the s-domain, and the equations that are given usually aren't explained how they came about. This makes it somewhat hard to understand the relationship between the math and circuits. As far as I understand, the sources generally have equations readily available for simple circuits, and rely mostly on theory of parallel and serial filters to derive larger circuits — is this correct? If so, is this the be all end all method of creating larger circuits, or what is the "correct" approach?

My goal with this is to 1) gain a better understanding of designing analog filters, and 2) take my own filters and make circuits for say guitar pedals and similar applications.

To be honest, I'd love to take any shortcut methods (minimal reading, maximum hands-on and results for motivation), such that I can get some circuits going right away, and then later take a deeper learning approach later on. Similarly, any "catches" or "gotya's" that I need to know about circuits, such that I don't end up blowing myself up, would be nice to know about ahead of time.

What is your recommended learning route to get started as fast as possible? Do you have any recommended material (articles, books, videos, etc.) for a speedy start?


texts: Why is my laser beam intensity profile not Gaussian?

texts:
I am kind of new to laser physics and signal processing, I hope I am in the right place to ask this question. I am also sorry for my not-great English.

I have been asked to study the transverse intensity profile of a laser beam produced by a $\rm He$-$\rm Ne$ source using a commercial webcam. I am placing the webcam sensor (without any lens) in front of the source and taking a snapshot. Then I am assigning each pixel an intensity value that is simply the recorded value in the red channel ($0-255$). If I choose a row of pixels and plot such value as a function of the pixel, I should expect a Gaussian profile. However it is not what I am observing.

In red you can see a Gaussian fit of the data. The actual data is shown in blue. I cannot really call this a good Gaussian shape, although it is some kind of a bell.

What could cause such a shape?
Could it be due to the laser source which could not be perfectly single-mode?
Or could it be due to the sensor?

It is interesting to note that I have also produced an image after making the beam pass through a single mode optical fiber, obtaining similar results; this makes me think that the problem lies in the sensor. However, the only modification the sensor causes and that I am aware of is Gamma correction, and as far as I know, that means that the red (or any channel) value is elevated to a power which is usually $<1$, but that should not change the actual shape. A Gaussian function elevated to a power is still a Gaussian function, although its parameters are different. Can anyone help?


texts: How do you "split apart" the output of a linear estimator?

texts:
Let's say I have a true signal $x$ which is corrupted by noise $n$ such that:
$y = x + n$
I am attempting to use a Weiner filter to estimate the true signal $x$. This estimator is of the form:
$\hat{x} = Gy$
Now let's say I keep the original equation and variables and add a known input bias $c$
$y + c = x + c + n$
The Weiner filter estimator would now take on the form
$\widehat{x+c} = G(y + c)$
The issue with this form is that I am now estimating $x + c$ instead of just $x$. Assuming I know the input bias value $c$, is it "okay" to simply break up the estimated output like the following?
$\widehat{x} = G(y + c) - c$
My confusion primarily comes about breaking up the output of an estimator. I assume because the Weiner filter is linear, this is okay.


texts:
I have been trying to understand the definition of the DCT and found a really good website describing its definition.

Here is the link DCT.

What I don't understand is the part of the DFT. Why is the summation of the DFT multiplied by $\frac{1}{\sqrt{2N}}$?

Thank you.


texts: How do Anti-Aliasing-Filters in audio signal processing work?

texts:
I am a little confused about the use of anti-aliasing filters. As on how I am confused, please consider the following task:

Consider a simple sampling rate conversion system with a conversion rate of 4/3. The system consists of two upsampling blocks, each by 2, and one downsampling block of 3.

The frequency spectrum has periodic repetitions at integer multiples of the sampling frequency. Therefore, upsampling creates additional - unwanted - spectral images. These can be cancelled by an anti-imaging filter:

The spectrum within the blue dotted rectangles should be my output Y_1. This output will be upsampled a second time, which I think results just having 4 spectrum repetitions until 2*pi.

But, how does downsampling work? I am given the following figure from my lecture notes:

So, apparently before downsampling we apply the anti-aliasing filter, which will result in the part of the spectrum within the blue-dotted shapes. Now, we increase Omega by a factor of 2, such that.. I don't know. What exactly happens here ?


texts: Unknown reason why the frequency spectrum of a "knocking force" (impulse) gets some spike

texts:
During my experiment (vibration analysis), I deliver a knocking force from a shaker (and the shaker itself has an impulse or an half-sine wave signal). The system automatically calculates with the FFT algorithm and gives me the following graph for the input force signal. Note that I was investigating different parameters, and this is not the proper experiment yet. Though the input signal would be probably kept the same for the main experiment.

My question is that: why do we spikes here, in the frequency spectrum?
Thanks for your help.


texts: change convolution function if signal is above threshold?

texts:
Let's say I have a signal $s(t)$ and two filters $f_1(t)$ and $f_2(t)$, I also have a threshold $A$.

Now I define $a_1(t)$ as $a_1(t) = \min(s(t),A)$, then I do the convolution of $a_1* f_1$.

Now I also define $a_2(t)$ as $a_2(t) = \max(0, s(t)-A)$, (ie. whatever is above $A$) then I do the convolution of $a_2*f_2$.

The total result $a$ is $a=a_1*f_1+a_2*f_2$.

Is there a name for such a kind of filter?

Basically, I do some convolution filter on some part of the input, and some other convolution on the remaining part of the input, then add the results

This can also by generalized with more than 2 filters, or continuous thresholds and continuously parametrized convolution functions


texts: How to compute the statistics of the DFT of correlated noise?

texts:
I would like to know how to compute the statistics of the discrete Fourier transform of a noise signal. To illustrate what I mean, I will first explain in detail a computation I have managed to do myself.

Suppose we have a discrete time series of values $x_n$ with $n$ from 0 to $N-1$. Each $x_n$ is a random variable, uncorrelated with the others, and Gaussian distributed with width $\sigma$. If I define the discrete Fourier transform

$$X_k = \frac{1}{N}\sum_{n=0}^{N-1} x_n e^{-2 \pi i n  k / N}$$

then I find that $X_k$ is a complex random variable with real and imaginary parts Gaussian distributed with width $\sigma/\sqrt{2 N}$. I did the computation by using the fact that the distribution of a sum is the convolution of the distributions, etc.

Now I want to know how to do this computation in the case that $x_n$ are correlated. How does one approach this problem? I can make the assumption that the process is Markovian.

I had originally asked this on the Computation Science site, but I think here is a better fit.


texts:
Is there a mathematical approach to initializing an IIR filter so that turn-on transients are minimized?
Background: I have a high pass Butterworth design that overloads at the start of a sine wave signal.  If I know the sine turns on at phase =0,i.e. at sine(t) = 0, can I preset the filter registers to minimize the overshoot?


texts: Basic difference between Fourier transform and laplace transform?

texts:
I have read few links about difference between Fourier transform and Laplace transform but still not satisfied

Please correct me if i am wrong
Simply put, the main difference between Fourier transform and Laplace transform is that real part is set to zero in Fourier transform while real part is non zero in laplace transform?


texts: What DSP techniques can be used to change the timbre of complex tone?

texts:
I'm experimenting with applying unusual 'spectral mangling' type special effects to audio.

1)  What DSP techniques can be used to change the timbre of a complex tone, such as that produced by the voice, or a musical instrument?

Are there any ways to manipulate the level of specific frequencies in the spectral profile of a complex tone. E.g. Assuming a tone with a fundamental frequency of 440 hz, use some DSP method to augment, for example the 2nd harmonic overtone at 1320Hz Khz, or remove another peak at say, 3Khz).
What's the best way to remove specific frequencies, or as narrow bands of frequencies as possible using DSP?

I use mixing/mastering, offline, and real time DSP VSTs, standalone programs, outboard effects units, but am interested in any methods which could work.

2) What alternatives (if any), are there to using techniques which introduce significant ringing or other artifacts which degrade the sound?

3) Are there any ways to 'clean up' audio after processing, removing artefacts which have been introduced?


texts: Beginner level : Help with terminologies : smooting length, length of channel, channel equalization delay, blind system identification, equalization

texts:
Say the channel model is an univariate FIR filter with true coefficients $h=[1,a_1, a_2 ]$. I am learning algorithms for system identification and channel equalization. For this, I am implementing programs provided in Matlab. One such Algorithm is Constant Modulus Algorithm for blind channel equalization.

Code is the link for the code : lines 80-82 present the variables $L,ChL,EqD$  which denotes the smoothing length, length of the channel & channel equalization delay respectively. I have confusion regarding the terms. Shall appreciate explanation with an example pertaining to my case.

1. Smoothing length : In text books, $L$ stands for the channel order, which for my model, Moving Average, MA(3), the order then = 3 ( I guess including the 1). FOr Autoregressive model AR(p) $p$ is the model order. I am confused if smoothing length has got anything to do with Kalman Smoother ? What is the difference, if any, between smoothing length and channel order? What will $L$ be for the FIR model that I have?

2. Length of the Channel, ChL = Is this related to number of data samples ?

3. Equalization delay, EqD = In general, I was under the impression that the model order is determined by the number of delays. So for MA(q) : $q$ is the number of delays, similarly for AR(p) : $p$ is the number of delays. Then what is Equalization delay?

When considering blind algorithm, we assume no knowledge of the process noise that drives the model. But, When we assume that the model order is known, then does it still become blind or non-blind estimation / equalization technique?

Please correct me where wrong.


texts: Continuous Wavelet Transform with Scipy.signal: what is parameter "widths" in cwt() function? How do time-frequency?

texts:
I search to display a time-frequency signal with an original discrete temporal signal (sampling step = 0.001sec). I use Python and the library Scipy.signal. I use the function cwt(data, wavelet, widths) to do a continuous wavelet transform, with the complex morlet wavelet (or gabor wavelet).

First step: Obtain a scale-translation signal.
In doubt, I associate directly the array “widths” with the array of the possible different scales. Because, I don’t understand what is parameter width if it’s not scale. Perhaps, you will tell me “it’s the width of your current wavelet”! But, even if it is, I don’t know how linked width with scale…

My second problem is to find and display the equivalent with frequency. In literature, I find this formula: Fa = Fc / (s*delta), where Fa is the final frequency, Fc the center frequency of a wavelet in Hz, s the scale and delta the sampling period. So, ok for scale (if I find the link with the width) and delta (=0.001sec), but it’s more complicated with center frequency of the wavelet. In scipy documentation, I find that:
“The fundamental frequency of this wavelet [morlet wavelet] in Hz is given by f = 2*s*w*r / M, where r is the sampling rate [s is here Scaling factor, windowed from -s*2*pi to +s*2*pi. Default is 1; w the width; and M the length of the wavelet].”
I think it’s the center frequency, is it? Maybe the solution to my first problem is here also (scaling factor and scale..?)..?

Thank you


texts: Understanding the attenuation of a signal

texts:
I have been told that the signal in the graph below was obtained by passing a seismogram through a high-pass filter to obtain the high frequency-small amplitude waves common in tremors far from the epicenter of an earthquake.

The high-pass filter used to unmask the tremor was selected to avoid the 'attenuation' of the tremor signal.  What exactly is meant by 'attenuation' in this context?


texts: finding the output of these three systems when connected in series?

texts:
I have just started to study Oppenheim's "Signals & Systems, Second Edition" and there is this easy-looking problem that has evaded me its solution for the past 48 hours:

consider systems $S1$, $S2$, and $S3$, where

$ x[n] \xrightarrow{S1} y_{(1)}[n] = \begin{cases}
x[n/2] & \text{if n is even;}\\
0 & \text{if n is odd.}\\
\end{cases}$

$ x[n] \xrightarrow{S2} y_{(2)}[n] = x[n] + (1/2)x[n-1] + (1/4)x[n-2] $

$ x[n] \xrightarrow{S3} y_{(3)}[n] = x[2n]$

To work out the result of $x[n] \xrightarrow{S1,S2,S3} y_{(1,2,3)}[n]$, I first tried to work out the result of $x[n] \xrightarrow{S1,S2} y_{(1,2)}[n]$ as follows:

$x[n] \xrightarrow{S1,S2} y_{(1,2)}[n] = \\
y_{(1)}[n] + (1/2)y_{(1)}[n-1] + (1/4)y_{(1)}[n-2] = \\ \begin{cases}
x[n/2] + (1/4)x[(n-2)/2] & \text{if n is even;}\\
(1/2)x[(n-1)/2] & \text {if n is odd.}\\
\end{cases}$

after that, feeding $y_{(1,2)}[n]$ to $S3$, we have:

$y_{(1,2)}[n] \xrightarrow{S3} y_{(1,2,3)}[n] = \\
y_{(1,2)}[2n] = \\ \begin{cases}
x[2n/2] + (1/4)x[(2n-2)/2] & \text{if 2n is even;}\\
(1/2)x[(2n-1)/2] & \text {if 2n is odd.}\\
\end{cases}$

which simplifies to $y_{(1,2,3)}[n] = x[n] + (1/4)x[n-1].$

but the solution manual says the output is in fact $x[n] \xrightarrow{S1,S2,S3} y_{(1,2,3)} = x[n] + (1/2)x[n-1] + (1/4)x[n-2].$ :(

What I want to know is, who is actually wrong: me or the manual? Thanks a lot in advance!


texts: How can I find the transform function, amplitude and frequency response of a digital filter in python?

texts:
I have applied a laplace filter mask to an image and now I want to find the amplitude and freqency response of a laplacian filter: [[1,1,1], [1,-8,1], [1,1,1] ]. I know I need to first find the transfer function, however, I am unable to do this as well programatically in python.

Here is the code I have currently:

from matplotlib.pyplot import imread
import matplotlib.pyplot as plt
from scipy.ndimage.filters import convolve
import numpy as np
import cv2 as cv

laplaceMask = np.array([ [1,1,1], [1,-8,1], [1,1,1] ])
img = imread("download.jpg")

channels = []
for channel in range(3):
    result = convolve(img[:,:,channel], laplaceMask)
    channels.append(result)

img = np.dstack((channels[0], channels[1], channels[2]))

plt.imshow(img)
plt.show()


texts: Initial conditions for the LTI systems described as a difference equations

texts:
Why do we need the initial conditions to be zero for the LTI systems described as a difference equations?

First question is why do we need it for linearity? I can't think of any example of the non linear system described as a difference equation with constant coefficients.
The definition of TI system is that it does not depend on particular time the input is applied. I can't understand how does this definition relate to the initial conditions? Those initial conditions will be the same whenever we apply the input so why they have to be zero?

Thanx


texts: How do I calculate the SNR of a noisy signal without having the noise?

texts: I wish to calculate the signal-to-noise ratio (SNR) of a speech signal in MATLAB. I only have the speech signal, and the speech-in-noise, but I don't have the noise file alone. How do I calculate the SNR of the speech in noise?


texts: FFT subset output yields computational efficiency

texts:
The traditional FFT gives $N=2^k$ outputs for $N=2^k$ inputs (generally complex.)  The time effort is $O(N \cdot \log(N))$.  If all I want is a compact subset of $M \ll N$ output frequency points, is there some method to "depopulate" the FFT late butterflies and achieve more efficiency than the depopulation burden adds?

The question of the FFT-like algorithm for fast DTFT computation? is close, but first misses the mark being really an off-frequency question, then the Goetzel transform efficiency threshold is far too low for my case.


texts:
Let's say we have 2 signals of frequencies 4 kHz and 8 kHz.

If we sample both signals with sampling frequency 10 kHz and pass the signal through a low pass filter with cut off frequency of 10 kHz, what frequencies will be present in the output signal?


texts: Find the stomata in a plant microscopy image

texts:
Here is a question for image processing experts.

I am working on a difficult computer vision problem. The task is to count the stomata (marked below) in DIC microscopy images.  These images are resistant to most superficial image processing techniques like morphological operations and edge detection. It is also different from other cell counting tasks.

I am  using OpenCV. My plan is to review potentially useful features for stomata discrimination.

Texture classifiers

DCT (Discrete cosine transform/frequency-domain analysis)
LBP (Local binary patterns)

HOG (Histogram of oriented gradients)
Robust feature detectors (I am skeptical)

Harris corners
SIFT, SURF, STAR, etc.

Haar cascade classifier/Viola-Jones features

And possibly design a novel feature descriptor. I am leaving out the selection of a classifier for now.

What have I missed?  How would you solve this? Solutions for similar object detection problems would be very helpful.

Sample images here.

After bandpass filter:

Canny edge detection is not promising. Some image areas are out of focus:


texts: Vector Quantization in Hidden Markov Models

texts:
In the field of Gesture Recognition and Speech Recognition (using HMM classifiers), do multiple HMM classifiers (one for each class) share common Codebook (i.e common set of symbols) in Vector Quantization?

If that is not the case, how does the test/classification stage work?

To which Codebook the test input vector belongs to (we need this to determine the symbol, in turn which is fed to the HMM for classification)?

Some background info:

One of the most earliest works done in Action Recognition used different set of symbols for each of the HMM classifier.
This journal in speech recognition by Rabiner suggest using common Codebook (with different HMM classifiers) for identifying 10 digits uttered by large number of users (in page 12),


texts: What is the best way to determine the process noise matrix $\mathbf Q$ of a Kalman filter?

texts:
It seems like most of the resources online suggest to determine the values of the process noise matrix $\mathbf Q$ through trial and error. However sometimes trial and error doesn't work, so I would like to know whether there is a way to approximate it e.g. via maths or whatever.

To clarify, I am talking about the matrix $\mathbf Q$ which can be found in the following formula:

$$\mathbf M_k = \mathbf \Phi_k + \mathbf P_{k-1} \mathbf\Phi^T_{k} + \mathbf Q_k$$ Where $\mathbf M_k$ represents the covariance matrix (before an update in this case). I know that
$$\mathbf Q = E\left[\mathbf w \mathbf w^T\right]$$

where $\mathbf w$ is a white-noise process. Which is part of this classical formula:

$$\dot{\mathbf x} = F\mathbf x + G\mathbf u +\mathbf w$$

But I don't know how to determine $\mathbf w$ neither...
To the people who are familiar with this very nice book (Fundamentals of Kalman Filtering, Zarchan et al.) you can find the former equation on page 131 and the second one on page 129.


texts: Using oversampling to increase resolution of a DC-signal as input

texts:
Currently I'm working on a project which uses oversampling to increase the resolution of a 12 bit ADC to a maximum of 16 bits. My goal is to fully understand the theory behind oversampling and why it is increasing the resolution. As far as I understood this topic, oversampling and decimation increases the resolution because the white noise of the input signal is distributed along a larger frequency span. After the oversampling the signal gets low pass filtered, so that we achieve fewer noise in the frequency span of interest (see the first figure).

So far so good, but I think i'm still not 100% confident of what is going on here:

How is the low pass filtering achieved? I know that I need to sample the input signal $4^n$ ($n$ = additional bits) times for every additional bit. After adding all samples together, the sum gets right shifted by $n$, which is equal to divide the sum by $2^n$. Is this the low pass filtering or do i need to somehow low pass filter it after i shifted the sum? In the application note AN118 that I used to understand this topic it seemed like the process of shifting is the low pass filter.

My input signal is a dc signal (sensor output) and I'm getting kind of a headache to understand why oversampling can increase the resolution of dc signals. Inside the appendix A of the application AN118 is a good approach to explain the reason for the increased resolution with oversampling. It also shows where the often used equation f_os = 4^n * f_s (f_os = oversampling frequency; f_s = sampling frequency) comes from. My problem is this whole appendix refers to AC-Signals... You can see this with for example equation 8 inside the application note, where the in-band noise power is calculated with the integral from 0 to f_m (f_m = highest spectral component of the input signal). In case of a dc signal as input i got f_m = 0. Am I missing something or can someone maybe explain to me why the resolution of a dc signal can be improved with oversampling?

Is something like SNR a good specification to describe the quality of a dc signal? I've read from a source that SNR is an AC specific attribute, for example in Table 2 of this Texas Instruments document. This document also states that ENOB can't be calculated for dc signals and that i need to calculate the "effective resolution" for dc signals. This was the first time i've read something about this... Can someone verify that for sure dc signals can't be specified with SNR and ENOB?

At the end of my project I would like to compare different resolutions achieved with oversampling. What would be a good specification to compare the measured results? Maybe the "effective resolution" or the overall noise calculated with the variance?
I hope someone can answere my questions. At the moment i'm struggling a bit to fully understand this topic...


texts:
$$h(t) = \frac{1}{1 + t^2}$$ and is it IIR or FIR filter. I tried finding the Laplace transform of this filter to get the data flow diagram with 5 taps and T=2s, however, I am unable to solve this. Probably the solution is simpler than I tried, however, being novice in DSP, I am unable to solve the same.

Thank You.


texts: How to use a FIR on a signal with a drifting central frequency?

texts:
I have been working on this subject for a while and I couldn't get a proper answer.

I am trying to estimate the phase of the input signal :

$x(t)=Ae^{j(at^2+bt+c)}+n(t)+i(t)$

Where $n(t)$ is a AWGN, $i(t)$ is an interfiring signal and $a,b,c$ are constant values.
As you can see, the central frequency of the useful signal changes over time ($f_c(t)=2at+b$).

First, I have to filter out the interfering signal so I implemented a low-pass FIR filter H and I compute the output of the filter as follow $y[n]=\sum{h[k].x[n-k]}$. However because of the moving frequency my signal will not stay in the filter passband very long.

One solution I had was to estimate the central frequency of the signal every 1 second and shift the input signal by this estimate (which means the I add to the FIFO of my FIR the samples $x[k]e^{-2\pi f_ck/F_e}$) but when I estimate the new central frequency I will have samples with different central frequencies in my FIFO and the filter doesn't like that (I can see phase jump at the time of frequency change). So can anybody tell me how I am supposed to tackle this issue ?

Thank you !

EDIT : Thank you very much for your answers, I realized while reading your answers that I simplified my problem too much. The actual evolution of the phase of my interest signal is much more complicated than a second order polynomial. My signal is more like
$x(t)=Ae^{jd(t)}+n(t)+i(t)$ with $d(t)$ a very complex function taking into account the distance and the speed of the object emitting the signal (which is a pure carrier). Also, the $i(t)$ is a function representing other object emitting other pure carrier at different speed and position (so my signal is more like $x(t)=\sum{A_ke^{jd_k(t)}}+n(t)$). I want to get the phase of each signal separately so first of all I did a small accumulation of signal to compute the spectrum of the signal and get the approximate frequency of each signal. Then I apply several FIR filters, one for each signal hoping I would get each $A_ke^{jd_k(t)}$ separately. The problem I have is that the frequency of each signal moves a lot and I have to "adapt" my FIR filter (and I don't know how to do that...)
I talked about a second order phase polynomial because it is a good model of the function $d(t)$ on a short term observation.

Thank you once again for the time and effort you put in your answers.


texts: Real-valued ringing when zero-padding odd-length FFT

texts:
So I'm trying to write a frequency-domain interpolator that zero-pads the frequency response of a signal and inverse transforms.  There's two cases I have to deal with:

Even-length response - have to split the $F_s/2$ bin because it's ambiguous.  So I copy the negative part of the spectrum, and add n*(interp-1)-1 zeros in between.
Odd-length response - there is no $F_s/2$ bin so just split positive/negative frequency and insert n*(interp-1) zeros between them.

The code that does the zero-padding can be seen here

// Copy negative frequency components to end of buffer and zero out middle
//  inp    - input buffer of complex floats
//    n    - transform size
//  interp - interpolation amount
void zero_pad_freq(cfloat_t *inp, size_t n, size_t interp) {
    if ((n % 2) == 0) {
        memmove(inp + n*interp - n/2, inp + n/2,     n/2*sizeof(cfloat_t));
        memset (inp + n/2 + 1, 0,       (n*(interp-1)-1)*sizeof(cfloat_t)); // Duplicate Fs/2 so we need one less zero

        inp[n/2]          /= 2.0;
        inp[n*interp-n/2] /= 2.0;
    } else {
        memmove(inp + n*interp - n/2, inp + (n+1)/2, n/2*sizeof(cfloat_t));
        memset (inp + (n+1)/2, 0,         (n*(interp-1))*sizeof(cfloat_t));
    }
}

The first case is working fine, I'm testing it on a chirp signal and it interpolates just fine, there's a little numeric noise, but it's round tripped through an FFT so what can you do (first $50 \mu s$ or so of the signal show):

The problem is with the odd-length transform, I'm getting a pretty heinous transient response on the real samples only ($50 \mu s$ again, real):

The imaginary channel has a small ripple on it, but not nearly as bad:

It's like I've screwed up my $F_s/2$ bin in the odd case, but there is no $F_s/2$ bin, so I'm very puzzled.  Anyone have any thoughts?


texts: Should we use all $N$ correspondence points to compute Homography or just $4$ points?

texts:
Given correspondences of $(x,x')$, we want to find a Homography matrix $H$ that maps $x' = Hx$. If we use 4 corresponding non - collinear points and the matrix $H$ has rank $n-1$, then the homogenous equation $Ah =0$ has a unique solution up to a constant. Suppose we find more than $4$ corresponding pairs of points, we get an overdetermined system of equation and we have to use least squares minimisation to find the solution.
I am curious as to whether it is better to just use $4$ points or $>4$ points to find the Homography ?
Also, it makes sense that from a linear algebra perspective, I am unable to solve the homogenous equation exactly when we have an overdetermined system. However, from a planar perspective, given my correspondence points are good, shouldn't there be an exact solution ?


texts: Is there any matlab (default) function for Discrete Hartley Transform? And is it applicable for audio signal? I want to apply it instead of DWT or FFT.


texts:
Is a passive LTI system will always be a stable(bibo and lyapunov) system?
Or in other words
Passivity of a LTI system implies stability?
If above statements are true then there must be some mathematical relation (description) should be there between passivity and stability, so anyone can help me about its mathematical description (if it exists)
And if passivity doesn't implies stability then there must be some counter example to it to prove my assumption wrong.
So whether my assumption is wrong or there exists a mathematical relation between passivity and stability?
I asked this question here
https://electronics.stackexchange.com/questions/513052/stability-of-passive-networks
but I didn't find any useful answer so I posted it here


texts:
I would like to remove vertical line from an image (an example). I took a 2D FFT and try to apply a mask to suppress the line. Nonethelesse the approch is not very efficient, because i lose an important part of information. How can i improve the treatment of FFT data? In FFT, how find the line?

My piece of code :

import numpy as np
import matplotlib.pyplot as plt
from skimage import io
from skimage import data, img_as_float

Path_input = "C:\\Users\\yoyo\\Desktop\\"

imggray = img_as_float(data.astronaut())[:,:,0]*255 #opening image
imggray[:,254:255] = 0 #force a vertical line
plt.imshow(imggray);plt.show()

imfft = np.fft.fft2(imggray)
mags = np.abs(np.fft.fftshift(imfft))
angles = np.angle(np.fft.fftshift(imfft))
visual = np.log(mags)
visual2 = (visual - visual.min()) / (visual.max() - visual.min())*255
plt.imshow(visual2);plt.show()

mask = io.imread(Path_input + 'mask_astro.png')[:,:,0]
mask = (mask < 100)
visual[mask] = np.mean(visual)

newmagsshift = np.exp(visual)
newffts = newmagsshift * np.exp(1j*angles)
newfft = np.fft.ifftshift(newffts)
imrev = np.fft.ifft2(newfft)
newim2 = np.abs(imrev).astype(np.float32)

plt.imshow(newim2);plt.show()

---- EDIT LATER ----

My "real" image with darkest parts that i have to maintain and stripes of noise.

and you can download here : https://image.noelshack.com/fichiers/2018/36/5/1536312281-test.png


texts: Mathematical Approach to Detect If a 2D Signal Is Separable

texts:
In the Special 2-D Sequences page the following examples are demonstrated,

2D dirac
2D diagonals
2D unit step function

Is there a more defined method or series of steps for determining if a function is separable or not? Other than:

$$x(n_1, n_2) = f(n_1)g(n_2) $$

Note: is there any difference between this and saying that:

$x(n_1, n_2) = f(n_1)\times g(n_2)$, where $\times$ is the orthogonal product (See SIMG-716 Linear Imaging Mathematics I section 2.1)? Is this simply referring to the inner product?

To me the solutions are highly intuitive and the requirement to solve such problems solely relies on the understanding of properties associated with the function.


texts: Suppose we have an unstable all-pole filter with transfer function $H(z) = A(z)^{-1}$. What is the best way to design a stable filter with frequency response as close to that of $H$ as possible. (I don't really know enough dsp to better define ''as close as possible'', but I am thinking in $L^2$ sense)


texts: Compression Sensing for Blind Source Separation

texts:
I am new to Signal Processing, and am interested in compression sensing for audio files. CS is based on the algorithm that, given some sampling of a signal $f$ in order to obtain a smaller (compressed) signal $b$ such that:
$$b = \phi f$$
such that
$$f = \Psi c$$
and $f$ can  be expressed by a linear combination of basis functions $\Psi$ and its coefficients $c$ by doing a nonlinear constrained optimization on the L1 norm of  such that :
$$\phi f = b$$
I understand how $f$ can be extracted from a sample without reaching the Nyquist limit, but I'm interested in the application of Compressed Sensing to blind source separation. With blind source seapration, to use standard notation
$$A = \phi\Psi , c = x  $$
Thus, L1 norm minization can be represented as
$$Ax = b$$
I have looked at Cleve Moler - Magic Reconstruction: Compressed Sensing for an example How can we adapt compressed sensing when we don't know A (since the whole point of BSS is that you do not know the mixing matrix) and we have to reconstruct more than one signal from b?


texts: How up sampling of a signal is done in real world?

texts: I know that up sampling by L means padding L-1 zeros between each samples. My question is how it is done practically? Do we use MATLAB or any other software, or there are hardware circuits to do so.


texts: Random process $X(t)$ with autocorrelation function given find the mean and the variance

texts:
Autocorrelation function is
$$R_{xx}(\tau)=\frac{20}{1+2\tau^2}$$
So at $\tau=0$$$R_{xx}(0)=20=E[X(t)X(t)]=E[X^2(t)]$$ The variance is
$$\mathrm{Var}[X(t)]=E[X^2(t)]-E^2[X(t)]=20-E^2[X(t)]$$ As $X(t)$ is WSS (wide sense stationary) the mean is a constant. Is there any way to find its numerical value?


texts: I want to plot a gaussian pulse in MATLAB in the frequency domain, with its centre at the highest value

texts:
I have tried one way but i am getting 2 peaks at the edges as shown in the picture. I have used 1 part of my code for the signal generation where i used the gauspuls command (reference: www.simsonic.fr) and for plotting the frequency spectrum, i picked up the code from a previous post.

The screenshot contains the relevant signals.the top for the frequency domain signal, bottom one is of the source signal.

%for signal generation

f0=0.5;
t0=1.5;
bndwdth=2;
duration=2*t0;

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

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

%for fft:
midf=5e5;
figure(1)
L=1000;
fs=midf*10;
%NFFT = 2^nextpow2(L);
NFFT=2286;
Mb = fft(signal,NFFT)/L;
f = fs*linspace(0,1,NFFT); %Frequency Vector
plot(f,abs(Mb))
title('Gaussian Pulse in Frequency Domain in bone')
xlabel('Frequency (f)')
ylabel('|Magnitude(f)|')


texts:
I am trying to figure out if there is anything more to the phrase 'complex band shifting' than just, a simple band shift. Disclaimer, this is NOT a communications problem, and I only say this so that some answers do get carried away with bit rates, etc.

Basically, let us say that I am interested in a signal, that exists in the passband, from frequencies $f_1$ through $f_2$. Let us further assume that I band pass filter this signal with a brick wall BPF. Then, I multiply it with a sine wave of frequency $f_1$. This operation shifts my signal 'down to baseband', and at base band, it now has a (single-sided) bandwidth of $f_2$ - $f_1$, (as it did before). (I can now LPF to get rid of the double frequency component). That is, my signal now exists from $0$ to $f_2$-$f_1$.

To me, this is all there is to band shifting. That is, the key component is a multiplication with a real sinusoid.

Why then, and what does, the phrase 'complex band shift' entail exactly?

I understand that a multiplication with a complex sinusoid is done in comm systems because we then try to decipher phase from the I and Q, so I get that, but if I simply want to shift bands of a signal, (even to an arbitrary frequency), what is 'complex band shifting' and is that always necessary? Cant I just use a real sinusoid for general band shifting?

Thanks.


texts:
i am working in digital image restoration field, recently i have studied image degradation model which is like that.
g(x,y) = h(x,y)*f(x,y) + n(x,y)

where

g(x,y) is the degraded image

f(x,y) is the image that we wish to capture or its is latent image

h(x,y) is image degrading function, and * is convolution operator

n(x,y) is noise

my question is that, can i change this particular model like that.

g(x,y) = h(x,y)*{f(x,y)+n(x,y)}


texts: Evaluate frequency reponse of cascaded filter with bitwidth truncation in matlab

texts:
How can we evaluate the frequency response of a series of fixed-point filters (H1,H2,...) cascading together in Matlab, such as:
12 bits input -> H1 -> bit truncation to 12 bits -> H2 -> bit truncation to 12 bits -> 12-bit output


texts:
I have a very simple question.

In Oppenheim book, it says that:

  If CT Fourier transform of $x(t)$ is $X(j\omega)$ then, CT Fourier
  transform of $x^*(t)$ is $X^*(-j\omega)$.

What I can't understand is what is $X^*(-j\omega)$? Or what is $X^*(j\omega)$? Can you tell me what I should understand when I see $X^*(-j\omega)$ or $X^*(j\omega)$?

For example $x(t) = (a + bj)t$ and $X(j\omega) = (c + dj)\omega$ what should be $X^*(-j\omega)$?

Thanks


texts: Pitch Detection / Harmonic Product Spectrum - What's wrong?

texts:
I've whipped together a small example that tries to use the Harmonic Product Spectrum Algorithm to extract the pitch of a simple sine wave. I'm not sure about my implementation and if i understood everything completely, that's why i'm hoping for some helping hand here.

The code is here: https://gist.github.com/akuehntopf/4da9bced2cb88cfa2d19
(You'll need JTransform and JavaFX to compile, though).

Please see what happens with the data during the processing using the
following graphs. Note: I'm using a sine wave with 160Hz and a sample rate of 8000. The duration of the sound data is 1 second.

First the sine wave is displayed (see first image).

Second, i'm applying a hamming window to the data (second image)

Then i calculate the power spectrum of the data as per

private float[] powerSpectrum(float[] window) {
    float[] powerSpectrum = new float[window.length];
    float[] fftBuffer = new float[window.length * 2 + 1];
    System.arraycopy(window, 0, fftBuffer, 0, window.length);
    FloatFFT_1D fft = new FloatFFT_1D(window.length);
    fft.realForward(fftBuffer);

    for (int i=0; i < fftBuffer.length / 2 - 1; i++) {
        float real = fftBuffer[2*i];
        float imag = fftBuffer[2*i+1];

        powerSpectrum[i] = (float)Math.sqrt(real*real + imag * imag);
    }

    return powerSpectrum;
}

Which gives me the spectrum in the third graph.

Observe that the peak is approximately at the right spot already, but we go for HPS, so continue (i'm not really sure about everything from here on):

Next, i'm multiplying the signal with its compressed form several times (with increasing compression)...

// 4. Compress
float[] spectrumCopy = new float[spectrum.length];
System.arraycopy(spectrum, 0, spectrumCopy, 0, spectrum.length);
for (int compression = 2; compression < 4; compression++) {
    for (int i = 1; i < spectrum.length; i++) {
        spectrum[i] = spectrum[i] * getCompressedSample(spectrumCopy, 1, compression, i);
    }
}

it uses the function

private float getCompressedSample(float[] buffer, int offset, int compression, int loc) {
    if (offset + loc * compression < buffer.length) {
        return buffer[offset + loc * compression];
    }
    return 0;
}

i found it somewhere on the interwebs. But my understanding is that we compress the power spectrum n times while taking only every 2, then only every 3 and so on samples from the power spectrum. The original spectrum is then multiplied with each of those compressed spectra.

The result i get is the fourth graph

Observe that the peaks have changed.

From my understanding the next step would be to find the bin with the highest peak, interpolate (i use quadratic interpolation) and recalculate the pitch frequency using this formula:

private float getFrequencyForIndex(int index, int size, int rate) {
    float freq = (float)index * (float)rate / (float)size;
    return freq;
}

However this gives me wrong results. I'm pretty much stuck. Any help with this topic would be very much appreciated! Thanks in advance!

For the next step i'm hoping to be able to extract pitch frequency data from a musical instrument (guitar/ukulele), but of course first things first :-)


texts: Mutual Information, is this calculation correct or should I contact the author about fixing it?

texts:
Background information

According to Cover's text book on information theory  Mutual Information is calculated as:

$$(1)  I(W;C) = \Sigma_k\Sigma_i P(C_k,W_i)log(\frac{P(C_k,W_i)}{P(C_k)P(W_i)})
$$

Now if we were interested in a single Word W_i we would just calculate this with one sum:

$$(2)   I(W_i;C) = \Sigma_k P(C_k,W_i)log(\frac{P(C_k,W_i)}{P(C_k)P(W_i)})
$$

OR we can rewrite this as

$$ (3)   I(W_i;C) = \Sigma_k P(C_k,W_i)log(\frac{P(W_i|C_k)}{P(W_i)})
$$

because
$$ (4)  P(W_i|C_k) = \frac{P(W_i,C_k)}{P(C_k)}$$ (unless i am incorrect)

Question
While reading a paper on feature selection using information theory measurements I came across A segment of the text:

Should the equation be more like (3), and I should contact the Author, or have I mis-understood/interpreted mutual information?


texts: How do I know if the given signal is discrete?

texts:
If I am given an expression of a signal then how do I know if it is a discrete signal or not?

If this signal is given how do I know if it's a discrete or continuous?

$$x(n) = e^{j (n \pi/2 + \pi/8) } $$


texts: How to perform FFT on stored digital data


texts: Apply decay envelope to waveform without multiplication

texts:
I'm trying to apply a decay envelope to a sinusoidal waveform using Verilog. Hardware constraints prevent use of multiplication to simply multiply by the envelope. The sine values and envelope values, for which I'm simply using $1-e^{-x}$, are stored in ROMs to speed everything up.

Right now I'm simply calculating the output sample by subtracting the exponential value from the sine value if the sine value is larger. If the exponential value is larger than the sine value, the output sample is zero. This works pretty well, but the current exponential value to subtract needs to be modulated by the current location on the sine wave so that the full value is subtracted when the sine is maximum and zero is subtracted when the sine is zero.

I'd appreciate any advice that anyone can give and I can provide more details if required.


texts: Why is level of power spectrum dependent on FFT resolution?

texts:
I created a sinusoidal wave in some noise, and plotted the power spectrum of the signal using two periodogram estimates (welch procedure). One estimate is 'high resolution' - i.e. it uses a longer window than the corresponding 'low resolution' estimate. The result and code used to generate it are shown below. For reference: the y-axis will have units of |signal|^2 and the x-axis are just FFT samples.

My question is: From the plot, it would appear that there is overall more power content in the low resolution (red) plot, although this is obviously not the case, the input signal is the same in both cases. This poses some questions about how we compare signals, it seems it is critical to use the same parameter values (like window lengths) when comparing the energy content of different signals.

For reference, I do understand that the length of the window in the welch estimate will determine the resolution of my spectrum, and that the number of FFT points used will determine the smoothness of my output, which has nothing to do with resolution.

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));

nWindow1 = 100;
noverlap1 = 0*nWindow1;
nfft1 = 1000;
[pxx1, f1] = pwelch(x,hann(nWindow1), noverlap1, nfft1, Fs, 'centered', 'power');

nWindow2 = 1000;
noverlap2 = 0*nWindow2;
nfft2 = 1000;
[pxx2, f2] = pwelch(x,hann(nWindow2), noverlap2, nfft2, Fs, 'centered', 'power');

figure
hold on
plot(f1, 10*log10(pxx1),'r')
plot(f2, 10*log10(pxx2),'b')
legend('Low Resolution Estimate', 'High Resolution Estimate')


texts: DFT, cepstrum, LPC for feature extraction

texts:
I am a newbie in speech processing and experimenting to get a feel. I have extracted some speech segments using a window function and I want to find distance between a pair of segments. First, I took DFT of a segment and downsampled from the resulting frequencies (like 16 out of 256) and used them as features.

Then, I saw cepstrum (cceps in matlab) and linear predictive coding (lpc in matlab) coefficients can be extracted from a speech segment. My question is; what is a good way to use these coefficients to find distance. Should I use the first k coefficients or downsample like in DFT? Any hints about their usage would be appreciated.


texts: Most natural interpolation to simulate increasing distance?

texts:
Background: I'm generating multiple datasets with varying scales of resolution to simulate taking pictures from varying distances and then comparing image classification performances. For example original images at 32x32, then downsampling to 16x16 and measuring the decrease in classification performance.

I'm using Matlab's imresize function to accomplish this. The docs say by default that it uses bicubic interpolation. Now I'm not sure if my following question even makes sense, but what I'm wondering is, is there a type of interpolation method that best mimics the decrease in resolution of an object that would happen naturally if you moved a camera farther away from said object? There's many different choices for the interpolation method parameter, and I'm wondering if I should choose a different one. Thanks!


texts: Multilevel partial wavelet reconstruction with PyWavelets

texts:
I'm looking for a way to partially reconstruct branches of a wavelet decomposition, such that the sum would recreate the original signal. This could be achieved in MATLAB using:

DATA = [0,1,2,3,4,5,6,7,8,9]
N_LEVELS = 2;
WAVELET_NAME = 'db4';
[C,L] = wavedec(DATA, N_LEVELS, WAVELET_NAME);
A2 = wrcoef('a', C, L, WAVELET_NAME, 2);
D2 = wrcoef('d', C, L, WAVELET_NAME, 2);
D1 = wrcoef('d', C, L, WAVELET_NAME, 1);
A2+D2+D1

ans =

    0.0000    1.0000    2.0000    3.0000    4.0000    5.0000    6.0000    7.0000    8.0000    9.0000

I'd like to achieve the same using pywt, but I'm not sure how to go about this.
The pywt.waverec function creates a full reconstruction, but doesn't have a level parameter for partial reconstruction.
The pywt.upcoef function does what I need for a single level but I'm not sure how to expand this for multiple levels:

>>> import pywt
>>> data = [1,2,3,4,5,6]
>>> (cA, cD) = pywt.dwt(data, 'db2', 'smooth')
>>> n = len(data)
>>> pywt.upcoef('a', cA, 'db2', take=n) + pywt.upcoef('d', cD, 'db2', take=n)
array([ 1.,  2.,  3.,  4.,  5.,  6.])


texts:
I am using RANSAC algorithm for homography estimation between pairs of images taken with cameras which do not have any translation between them (pure rotation and change of scale/zoom). It works well in half of the cases. The correct output looks like this:

The red lines are filtered correspondences and the quadrilaterals illustrate how the homography distorts perspective.

Sometimes, however, many bad cases happen, like these:

I already have a simple test in the RANSAC loop. It makes a simple quadrilateral (a unit square) and transforms it with sample transform. Then it looks whether the transform kept its convexity.

Still, however, bunches of concave quadrilaterals come out.

Do you have any idea how to properly test the homography, if it behaves "nicely" and filter out the incorrect solutions?

I found some code where they test that none of the three transformed points are colinear. But this does not seem sufficient as it will not filter out deltoids and other "invalid" quadrilaterals...


texts: transmitting data using phase modulation of two sine waves passed through distortion

texts:
I have two sine signals with frequencies f1 and f2, and phase p1 and p2.

I modulate the phase of these signals to transmit two symbols simultaneously, chosen from the set [0,1,2,3]. For example, I represent 0 by 45 degrees phase shift, 1 by 90 degrees phase shift, etc.

The signals will pass through an optical communication channel, potentially getting distorted. This means there could be intermodulation products and harmonics.

I plan on using an FFT to determine magnitude and phase changes of the signals. I am hoping to retrieve the symbols from the phase changes I detect using the FFT. I want the magnitude to stay constant while the phase is shifting to transmit the symbols.

The sampling rate I'm working with is 7800 hz. I will potentially have to fit 16 transmission channels (frequencies) below 1500 hz. At any time, two frequencies randomly chosen from the set of 16 frequencies may be present, transmitting symbols as explained above.

Is this going to work, and what are the key issues to keep and mind and to work around? Any practical pointers or references are highly appreciated as I'm a novice in phase modulation and encoding/decoding of information in this fashion.

Thanks for helping!
B


texts:
I am really struggling to see how the lecturer took this transfer function and produced the bode plot that I have in my notebook. It is not the plotting so much that is confusing, I just don't understand how he manipulated the transfer function into this form that he did and why it needed to be in that form in the first place when it looked solvable in its original form.

  $\textrm{(b)}\quad T(s) = 180\displaystyle\frac{s(1+0.01s)}{(1+0.05s)(1+0.001s)}$

Looking at the transfer function shown I would have thought that all one would simply need to do is evaluate the gain at $20\log(180)$ and then plot the bode plot using the two zeros and two poles.

Supposedly to plot the asymptotic bode plot for this function correctly one needs to get it into the following form:

$$3600\times\frac{\left(\frac s{20}\right)\left(1+\frac s{100}\right)}{\left(1+\frac s{20}\right)\left(1+\frac s{1000}\right)}$$

Sorry for the bad format but basically the answer provided has it rising to its pole at $20$ on the $71\textrm{ dB}$ line and going flat until it hits $100 \textrm{ rad/s}$ and then rising up again at $20\textrm{ dB}$ per decade until it meets its pole at $1000$. I don't understand what is going on and everything I have written hitherto is just my observation of how the lecturer demonstrated.

If anyone with any knowledge on the mechanics of plotting bode plots by hand could shed some light on this for me I'd greatly appreciate it.
Simon.


texts: MLE parameter estimation -- confusion regarding some terms in the pdf of complex normal r.v (Part 2)

texts:
This question is based on the application of the pdf which was an earlier question of mine asked here Confusion regarding pdf of circularly symmetric complex gaussian rv

If $v \sim CN(0,2\sigma^2_v)$ is a circularly complex Gaussian random variable which acts as the measurement noise in this model
$$y_n = A + v_n \tag{1} $$ where $y$ is the observation and $A$ is a scalar unknown value which needs to be estimated. I am having a slight confusion whether there will be a 2 in the denominator of Eq(3) and Eq(4) with the $\exp(.)$ term. Based on the answer in the link, there should be no sqrt term with $\pi$ in the denominator, if $v \sim CN(0,2\sigma^2_v)$. If $v \sim N(0,\sigma^2_v)$ then there is a sqrt term.

Can somebody please check if I have correctly written  out the log-likelihood? I think I am missing a 2 in the denominator of $\exp[.]$ term in Eq(3) but I am not quite sure.

Thank you for your time and help.

$$P_y(y_1,y_2,...,y_N) = \prod_{n=1}^N\frac{1}{2\pi \sigma^2_v} \exp \bigg(\frac{-{({y_n-A})}^H ({y_n-A})}{2\sigma^2_v} \bigg) \tag{2}$$

taking log
$$\ell   = -N\ln(2\pi\sigma^2_v) -  \frac{1}{\sigma^2_v} {\bigg[{[\sum_{n=1}^{N} {(y_n - A)}{(y_n - A)}^{\mathsf{H}}  ]}\bigg]}. \tag{3}$$
$$ = -N\ln(2\pi\sigma^2_v)- \frac{1}{2 \sigma^2_v}{\bigg[ \sum_{n=1}^{N}y_n y_n^\mathsf{H}
- 2 \sum_{n=1}^N y_n A\bigg]}
  - \frac{1}{2 \sigma^2_v}{\bigg[  \sum_{n=1}^N   {AA}^\mathsf{H} ] \bigg]} \tag{4}$$


texts: Power spectral density of $\left(x(t)\right)^2$?

texts:
The relation between  $x(t)$ and output $y(t)$ of a non-linear device is expressed as
  $$y(t) = (x(t))^2$$
  Let $x(t)$ be zero-mean stationary Gaussian random process with auto-correlation
  $$R_x(\tau) = e^{-a\lvert \tau \rvert} , \quad a>0$$
  Find the output PSD?

I tried using the formula of output $\textrm{PSD}_y = \lvert H(f)\rvert^2 \cdot \textrm{PSD}_x$ , but stuck in process of finding $H(f)$ of a non linear system squarer. Not able to think of any other idea.


texts:
Given a system which is either non-linear and/or time variant, Is there any systematic way to determine if this system is invertible?
For example how to determine invertibility of this system?
$$\int_{-\infty}^\infty x(5\tau)d\tau$$


texts: What is the multiplication rate in FIR filter

texts:
Consider a system implementing a rational sampling rate change by 5/7: for this, we cascade upsampler by 5, a lowpass filter with cutoff frequency pi/7 and a downsampler by 7. The lowpass filter is a 99-tap FIR.
Assume that the input works at a rate of 1000 samples per second. What is the number of multiplications per second required by the system? Assume that multiplications by zero do not count and round the number of operations to the nearest integer.
Thanks in advance!