texts: Implementing an L-System conversion from R2 to R1 for symmetric, square matrices

texts:
I have been reading for some time about L-Systems, and specifically the Hilbert Space filling curve.  I am interested in writing a function to convert upper-triangular matrix coordinates into an ordered, 1-dimensional space filling curve.

The usual Hilbert space filling curve is derived with a single, first-order curve:

Then

I have successfully implemented the mapping algorithm listed on Wikipedia,

//rotate/flip a quadrant appropriately
void rot(int n, int *x, int *y, int rx, int ry) {
    if (ry == 0) {
        if (rx == 1) {
            *x = n-1 - *x;
            *y = n-1 - *y;
        }

        //Swap x and y
        int t  = *x;
        *x = *y;
        *y = t;
    }
}

//convert (x,y) to d
int xy2d (int n, int x, int y) {
    int rx, ry, s, d=0;
    for (s=n/2; s>0; s/=2) {
        rx = (x & s) > 0;
        ry = (y & s) > 0;
        d += s * s * ((3 * rx) ^ ry);
        rot(n, &x, &y, rx, ry);
    }
    return d;
}

//convert d to (x,y)
void d2xy(int n, int d, int *x, int *y) {
    int rx, ry, s, t=d;
    *x = *y = 0;
    for (s=1; s<n; s*=2) {
        rx = 1 & (t/2);
        ry = 1 & (t ^ rx);
        rot(s, x, y, rx, ry);
        *x += s * rx;
        *y += s * ry;
        t /= 4;
    }
}

However, I need to add two more first order curves to this process: one that is a right angle and the other which is a wedge.  These two new curves do not need to be rotated or flipped, and in combination with the usual Hilbert 1st order curve, one can fill an upper triangular matrix (see my sketch below):

That is, the upper left of the red right-angle again becomes a right angle, the upper right becomes a Hilbert curve, and the lower right part of the right-angle becomes a wedge, and so on.

I would like to implement the above proposed L-System conversion, but every single tutorial on L-Systems is about STARTING AT ZERO and then DRAWING A LINE that fills the space without crossing itself (i.e. writing a string of steps to trace the whole space, instead of converting a subset of coordinates in the space to there 1-D counterpart).

Can anyone offer any intuition on how I might incorporate these new restriction into the existing code?  Are there any resources that explain L-System conversions (as opposed to drawing systems), and how to construct them in code??


texts: Transform coding gain of two channel Haar Wavelet transform (first stage)

texts:
Question

My interpretation

Work

I'm attempting to do this in MATLAB.

I know that for discrete Random Processes going into an LTI system the autocorrelation function of the output is $R_y[n] = h[n]\star h[-n]\star R_x[n]$ where $\star$ is convolution.
The two tap Haar filters are: h1 = [0.707 0.707] and  h2 = [0.707 -0.707] which are orthonormal to each other.

I also know that transform coding gain is defined this way:
$$G = \frac{1}{N}\frac{\displaystyle\sum_{i=0}^{N-1} \sigma_i^2}{\left(\displaystyle\prod_{i=0}^{N-1} \sigma_i^2\right)^{\frac 1N}}$$

where $N$ is the number of samples (3 in this case?) and $\sigma_i$ is the variance of the $i$th transform coefficient.

This is my code:

h1 = 1/sqrt(2).*[1 1];
h2 = 1/sqrt(2).*[1 -1];
y1 = conv(x,conv(h1,fliplr(h1)))
y1 = 4.6500   14.3000   19.3000   14.3000    4.6500
y2 = conv(x,conv(h2,fliplr(h2)))
y2 = -4.6500    4.3000    0.7000    4.3000   -4.6500

After decimation:

Y1 = 4.6500   19.3000    4.6500
Y2 = -4.6500    0.7000   -4.6500

Calculated coding gains:

GTC1 = 9.533/(417.3134)^(1/3) = 1.2757
GTC2 = -2.8667/(15.1358)^(1/3) = -1.1589

My questions:

Am I calculating $y_1$ and $y_2$ correctly? As I understand it these are supposed to be the autocorrelation functions for a WSS Gaussian random process (because WSS into LTI gives WSS out). After decimation, can i use $Y_1$ and $Y_2$ to calculate coding gain? I'm getting a negative coding gain for the second channel ($Y_2$), does that make sense?

EDIT: Just realized that the Haar filters need not be 2 tap filters. Now I am more lost than before.


texts:
Using peak search, I'm able to put the cursor on any of the several peaks on the spectrum analyzer display.  I am trying to do something similar in software, with the output of the FFT of the radio spectrum.  Using normal peak detect functions (such as those included in Scipy) does not seem to work.  I think that there must be a particular algorithm used for radio, where particular circumstances exist (relatively flat noise floor, etc.) , what is it?

EDIT:

Here is a screenshot:

The orange markers are where the peaks are thought to be.

How I made this was filtering and using scipy.signal.argrelextrema, though there are a few different techniques (some on stack overflow) that seem to be just as tricky.


texts:
I want to do noise shaping in a 100kHz, 16 bit application, so as to shift all quantization noise to the 25khz-50kHz band, with minimal noise in the DC-25kHz band.

I set up mathematica to create a 31-sample error filter kernel via reinforcement learning which works well: After a little learning, I can get about ~16dB enhancement of high frequency noise for about the same amount of reduction in the low frequency band (the central line is the unshaped dither noise level). This is in line with the "Gerzon-Craven" noise shaping theorem.

Now to my problem:

I cannot manage to shape the noise even more even after extensive learning, although the Gerzon-Craven theorem doesnt forbid it. For example, it should be possible to achieve 40 dB reduction in the low band and 40 dB enhancement in the high band.

So is there another fundamental limit I am running into ?

I tried looking into the Shannon noise/sampling/information theorems but after fiddling with that a while, I only managed to derive a single limit from it: the Gerzon-Craven theorem, which seems to be a direct outcome of the Shannon theorem.

Any help is appreciated.

EDIT: more info

First off, the filter kernel that produces the above noise shaping, Note that the most recent sample is on the right side. Numeric values of the BarChart rounded to .01: {-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29} (Not exactly the bar char but produces similar curve)

Another note about error feedback implementation:

I tried two different implementations of the error feedback. First I compared the rounded output sample to the desired value and used this deviation as error. Second I compared the rounded output sample to the (input+error feedback). Although both methods produce quite different kernels, both seem to level off at about the same noise shaping intensity. The data posted here uses the second implementation.

Here is the code that is used to calculate the digitized wave samples. step is the stepsize for rounding. wave is the undigitized waveform (typically just zeros when no signal is applied).

TestWave[kernel_?VectorQ] :=
 Module[{k = kernel, nf, dith, signals, twave, deltas},
  nf = Length@k;
  dith = RandomVariate[TriangularDistribution[{-1, 1}*step], l];
  signals = deltas = Table[0, {l}];
  twave = wave;
  Do[
   twave[[i]] -= k.PadLeft[deltas[[;; i - 1]], nf];
   signals[[i]] = Round[twave[[i]] + dith[[i]], step];
   deltas[[i]] = signals[[i]] - twave[[i]];
   , {i, l}];
  signals
  ]

The reinforcement method:

The "score" is computed by looking at the noise power spectrum. The goal is to minimize noise power in the band DC-25kHz. I am not penalizing noise in the high frequency band, so arbitrarily high noise there would not diminish score. I am introducing noise in the kernel weights to learn. Maybe, therefore, I am in a (very broad and deep) local minimum, but I consider this extremely unlikely.

Comparison to standard filter design:

Mathematica allows to generate filters iteratively. These can have much better contrast than 36 dB when their frequency response is plotted; up to 80-100 dB. Numeric values: {0.024, -0.061, -0.048, 0.38, -0.36, -0.808, 2.09, -0.331, -4.796, 6.142, 3.918, -17.773, 11.245, 30.613, -87.072, 113.676, -87.072, 30.613, 11.245, -17.773, 3.918, 6.142, -4.796, -0.331, 2.09, -0.808, -0.36, 0.38, -0.048, -0.061, 0.024}

However, when applying those in the actual noise shaping, they are (a) clamped to the same ~40dB contrast, (b) perform worse than the learned filter doing actually no noise attenuation.


texts: How to calculate bandwidth of a wire, "air" or other channels?

texts: I'm really curious of how engineers identify what is the frequency range that could be sent through "wire" as channel? Do they put wire in front of them and make test that consists of sending sinusoidal signals through it by gradually increasing frequency, till maybe wire suddenly "burns off" and they measure that as highest frequency of a "wire's bandwidth"? How do they test "air" as channel where EM waves travel?


texts: Why does the griffin lim algorithm work?

texts:
the griffin lim algorithm works to reconstruct audio from a magnitude spectrogram, but i don't understand how reusing the same spectrogram every iteration with noise can produce decent sound

How from the random phase is ordered to be heard well?


texts: How to calibrate a signal with respect to frequency

texts:
I have recordings from accelerometers (non-flat response with respect to frequency) and I'd like to convert my voltage time series into acceleration values.  Is there a standard method for doing this?

My proposal to achieve this is:

Break the signal into equally spaced segments overlapping with 50%.
Apply a fft (with appropriate window) to each segment and apply the calibration values to each frequency.
Convert the segment back into a time series by applying an ifft and reversing the window transform.  We not have calibrated segments of the time series, although they are likely disjoint between window segments.
Apply a triangle window to all overlapping segments and add them into a single time series.  The triangle window should result in smooth transitions between window segments, resulting in a continuous, calibrated time series.

Are there any problems with this approach?  The window choice is dependent on use case.  Are there other approaches to achieve the same goal?


texts: Gaussian Mixture Model-UBM based for image recognition

texts:
Scenario:

Open set recognition: Unknown individuals may appear on scene (CCTV).
Gallery = 10 individuals with a mug-shot image per individual.

Approach:
Gaussian Mixture Model-UBM based

Use video sequences of 20 individuals of random people that might not appear in the gallery or in the scene as UBM (universal background model).

ubm = gmm_em(datalist, numberOfMixtures, EMiterations, downSamplingfactor, parallelWorker)

'datalist' contain feature vector of each individuals that assigned as UBM (20xnumberOfvideoSequences) cells. Each cell contains [numberOfMixtures x numberOfFrames] double.

Adapt the ubm into each individuals (in gallery) to create GMM speaker model.

for I = 1:numberOfIndividuals gmm{s} = mapAdapt(spkFiles(I,:),ubm, map_tau, config) end

spkFiles contain feature vector of the mugshot images for each individuals in the gallery. (10 x 1) cells, which each cell contain [numberOfMixtures x 1] double.

The approach above is based on MSR toolkit for speaker verification. However, I'm not so sure about implementing it for image verification. Is it what I had described above is correct?
The way to prepare the datalist and spkFiles?
How to decide the number of mixtures? Can I use any feature extraction method such as PCA?


texts: Sub-Band Extraction Techniques for Complex (I/Q) Signal

texts:
There is an application which processes Complex I/Q signal of 100 MHz bandwidth and it is required to extract independent sub-bands from this bandwidth at different center frequencies and bandwidths. I wanted to ask what are the techniques which can be used to extract this sub-band data?
Edit: FFT can also be considered as a filter bank. Can the FFT information be used to create independent sub-band data using IFFT etc.?


texts: Is convolution distributive over multiplication?

texts:
Is there any formula or expansion for

$$a(t)*[b(t) \cdot c(t)]$$

$$a(t) \cdot[b(t)*c(t)]$$

where $*$ denotes the convolution?

By expansion I mean something like $a(t)\cdot[b(t)+c(t)]=a(t)b(t)+a(t)c(t)$.


texts: Removing parallel edges from binary image

texts:
This is a binary image. I want to remove the portion in the red region, and keep the one in the green region.

The red region has the property of double parallel edges, whereas green does not. If someone could suggest an algorithm it would be helpful.


texts: Ideal Reconstruction of Upsampled Signal

texts:
Problem: The signal $cos(2\pi14100t)$ is sampled at $F_s = 400 Hz$. It is then upsampled with a factor 3 and then reconstructed ideally with a new frequency $F = 500 Hz$.

I now want to find the new signal that is created.

Incomplete solution: With the frequency before sampling $\pm 141000 Hz$. After sampling the normalized frequencies:

$f = \pm {14100 \above 1 pt 400} \pm k = \pm (35 + {1 \above 1 pt 4})\pm k= \pm {1 \above 1 pt 4} \pm k $ per sample

Interpolation with factor 3 yields the normalized frequencies:

$f = {\pm {1 \above 1 pt 4} \pm k \above 1 pt 3} = $ 3 frequencies $\pm k?$

I know how to solve it past that to get the final $F= \pm f * F_s$. I.e. $y(t) = cos(2\pi f_1t) + cos(2\pi f_2t) + cos(2\pi f_3t)$.

So my question is how do I calculate the 3 frequencies? Would greatly appreciate a simple solution and not just an answer so that I actually understand this.


texts: How do I interpret these autocorrelograms?

texts:
Hello signal processing wizards.  I am young and feeble, and would like some help interpreting my spike-train autocorrelograms of six isolated neurons.  These histograms were created using the code pasted below (hist(k)).  Rather than creating an autocorrelation of the actual spike trains, this code correlates interspike intervals (see the following post for another example:

https://stats.stackexchange.com/questions/148865/how-to-interpret-this-autocorrelogram-graph

I am glad to have gotten the code working, but now I am struggling to interpret what the histogram actually means.  I am aware that the Y-axis shows how strongly the signal is correlated with itself at a particular lag (the X-axis).  However, I am not certain how I can relate these findings to the firing pattern of the neuron.  For instance, if I see a peak at say, ~ 750 ms, does this mean that with a lag time of 750 ms, the firing pattern is more constant/predictable, or does it simply mean that I have more inter-spike intervals of 750 ms?  If someone could clear up my confusion, that would be much appreciated.  On a related note, if anyone has any tips on how to choose an appropriate bin size and window, that would be awesome.

P.S. The X axis ranges from -1000 ms to 1000 ms.  Bin size is currently set to 50 ms.

time_of_spike= sparse([],[],[],1,tmax/dt);

clus_ms_Mult_10 = round((clusterTimes*10));%cluster time-stamps were in ms

while index < length(clusterTimes);
    time_of_spike(clus_ms_Mult_10(index))=1;
    index=index+1;
end

% % %finding the interspike intervals matrix
I=find(time_of_spike)/10;
i=size(I,2);
for m=1:i
    for n=1:i
        ISI(m,n)=(I(m)-I(n));
    end
end

%calculating the autocorrelation
b=1;
clear k;
for i= MIN:1:MAX
    k(b)=size(find(ISI(:,:)>=i & ISI(:,:)<=i+50),1);%if i want to change "bin
                                                    %size", I can change i+1 to i+largerInt
    b=b+1;
end


texts: How to calculate Noise and distortion of an image?

texts:
I need to make an application in iphone which needs to calculate noise, geometric deformation other distortions in an image. How to do this? I have done some image processing stuff with opencv + iphone. But I dont know how to calculate these parameters.

1) How to calculate noise in an image?

2) What is geometric deformation and how to calculate geometric deformation of an image?

3) Is geometric deformation and distortion are same parameters in terms of image filter? or any other distortions available to calculate an image is good quality or not?


texts: I do not know how to create sturct variable in Maltab from column vectors. I have three column vectors, Voltage1 (32768x1 double), Voltage2 (32768x1 double) and Voltage3 (32768x1 double). I would like to create 1x1 struct variable with 3 fields (Voltage1-32768x1 double, Voltage2-32768x1 double and Voltage3-32768x1 double). Please help.


texts: How to integrate and differentiate sampled signals?

texts:
I am working on a control system which is defined by the equation.

As you can see,the equation is in continuous state.

$T_i$ = integral time and $T_d$ = differential time which are both constants and I can figure them out (Not of interest for now).

In the real life control system that I am trying to develop, I have signals sampled at regular intervals of time.

Consider them to be x1, x2, x3, …

e(t) = set_point - xt

xt is measured value at time t.

set_point is the desired value.

I somehow have to implement the above equation as an algorithm i.e apply the integration and differentiation in the constants (signal values that are sampled at regular interval of time ). Some references that I saw mentioned that integration and differentiation are mere addition and subtraction of the sample values over time.


texts: What's the difference between the Gabor-Morlet wavelet transform and the constant-Q transform?

texts:
At a glance, the constant-Q fourier transform and the complex Gabor-Morlet wavelet transform seem the same.  Both are time-frequency representations, based on constant-Q filters, windowed sinusoids, etc.  But maybe there's a difference that  I'm missing?

Constant-Q Transform Toolbox for Music Processing says:

  CQT refers to a time-frequency representation where the frequency bins are geometrically spaced and the Q-factors (ratios of the center frequencies to bandwidths) of all bins are equal.

Time-scale analysis says:

  That is, computing the CWT of a signal using the Morlet wavelet is the same as passing the signal through a series of bandpass filters centered at $f = \frac{5/2\pi}{a}$ with constant Q of $5/2\pi$.


texts: How to remove frequency offset with RRC as pulse shaping filter

texts:
I am working on a project with QPSK modulation format and a RRC filtering at the transmitter (TX) side. The channel introduces a constant dropper shift ($f_d$) to the received signal. To make the channel ISI free, the same RRC filter is used on the receiver (RX) side.
To remove this $f_d$, a practical ML frequency offset estimator is used [1].
This method is derived by assuming this $f_d$ is small enough. Therefore, after some mathematical manipulation, the decision statistics can be narrowed down to
$$\max(\lvert \text{fft}(c(k)^**y(k)\rvert),$$
where y(k) is the received r(k) filtered by a RRC matched to the RRC at TX, and c(k) is the transmitted data symbol (sample). To eliminate c(k) in the statistics, the $(c(k)^**y(k) )^4$ is used before fft operation is performed. The upsampling rate is 10. I couldn't see any issues in this process. Besides, I also tried to do $((c(k)*y(k))^4 )^{\frac 14}$ in order to eliminate the effect on the magnitude of the filtered received signal, but this gives me a huge DC component.
Finally, the algorithm works very well without RRC filtering at TX and RX. It also works very well without this $f_d$. So the problem should be something related this pulse shaping part.
Could anybody see any problems in this process?
[1]. U. Mengali et al: "Synchronization Techniques for Digital Receivers".


texts: In Plain English, What Are the Within and Between Class Variances in Otsu's Thresholding?

texts:
For the past week I've been through numerous research papers and tutorials on Otsu's thresholding. Even wrote to the authors, but every single one of the explanations are in terms of "Class" and "variance". I'm not very good at statistics, so it's hard for me to picture what these classes and variances mean in terms of the actual image. I have a few questions (marked in bold) that lead to the main questions.

What I've understood until now:
1. We need to separate the image into two parts. The dark pixels and the light pixels. These are what are called the two classes. For this explanation could we please just call them "dark" and "light" instead of "class".
2. A histogram array of 0 to 255 positions is created, where each position holds the total number of pixels of intensity 0, 1, 2,...255. Let's say for intensities 0 to 9, the number of pixels are 20,5,23,6,3,23,5,8,67.
3. Probability for any of these intensities (assume intensity 1) are calculated as = 1 / (180 * 200), where 180 is the length of the image and 200 is the width.
4. The sum of probabilities of the dark and light pixels are the weights of the dark and light pixels.
5. Question: When we calculate the mean of the pixels, what exactly does it imply? Is it like taking an average of the probabilities? If it is like an average, then if the threshold is at position 100 of the histogram, why can't we simply calculate it as: (sum of intensities from 0 to 99) / (99) ? What extra benefit do we have by doing a complex calculation of mean?
6. Then we have variances of distribution of pixels from the mean, which I understand, but Question: when we take the sum of variances of the light pixels and sum of variances of dark pixels, what does it imply? Does it speak of the amount of variation of intensity present in the dark pixels vs the light pixels?
7. Question: What is a "within class variance" and a "between class variance" in plain English? I can see the formula, but it makes no sense. Is the within class variance just a way of saying "how much the dark pixels vary with respect to mid-darkest pixel plus how much the light pixels vary with respect to the mid-lightest pixel"?
8. Question: What is a "between class variance" in plain English? The formula itself is odd, because there is no variance, but just a sum of the wights and means. But still, what does it mean in plain English, in terms of the pixels?

ps: Kindly refrain from statistical/mathematical terms and please try to explain it in plain English the way I've put some explanations in quotes above. I (and many people in this world) have a hard time understanding statistical jargon because it was taught in a way that never allowed us to fully understand it.


texts: Aliasing correspondance with respect to DFT

texts:
Given a sampling frequency Fs lets say we plot the magnitude of the fft of a temporal signal $x$, for different frequencies above the Nyquist frequency, to show the effect of aliasing:

Fs=10;
t=0:1/Fs:1-(1/Fs);
f=6;%7,8,9, ... to see aliasing starting a 6 since 6 is just above the Nyquist freq
x=sin(2*pi*f*t);
freq_range=Fs*((0:length(x)-1)/length(x));
plot(freq_range,abs(fft(x)));

Can someone help to formulate mathematically to which lower frequencies a given freqency f > Fs/2 will be aliased to?

Experimentally, we can see that (for integer frequencies, dont know if it really makes sense to use fractionnal ones?):

6 ->  4

7 ->  3

8 ->  2

9 ->  1

10 -> 0 <- should be zero? but the spectrum goes crazy and there are peaks everywhere.

Why? Then one would think that this would be periodic i.e. f=11 -> f_aliased=1, etc. back to zero in decreasing order, then again some strange things happens at 15 then it goes on cyclically etc. What is the relation?

Also, w.r.t. spectral resolution ($\Delta f = F_s/N$, let $N$ be the nb of samples, i.e. length(x)), one could ask if the next frequency abouve the Nyquist is in fact $f_{Nyq}+\Delta f = (F_s/2)+\Delta f$ ? But in this example the spectral resolution $\Delta f$ = 1 right?


texts: Vibration signal analysis of non-stationary movement

texts:
my task is to develop a condition monitoring system of a machine module. For that I am using the vibration sensors placed on the machine housing. I want to focus on only one module of the machine, a rack and pinion system, in which the rack carries out one reciprocating movement (back and forward) per machine cycle. I have carried out some test measurements to see if I see the vibrations corresponding to the cycle of the module. I could not see any clear vibrations that I could trace to rack and pinion module. Here is the qualitative representation (The signal above is the vibration, the curve below is the rotational speed of the rack):

The rotational speed curve is just an assumption, since I do not have any speed sensor in place to measure it.

There are many rotating gears in the transmission of the machine and I need to extract only the vibration signal of my module. What I want to do is to install a proximity sensor next to the pinion teeth to be able to measure its rotation. I would then know when exactly does my rotation start and end, and I would be able to differentiate the signal flanks to get the instantaneous rotational speed of the pinion. Once I have this information I would be able to extract the parts of the signal when the cycle occurs and analyse them using STFT or Order Tracking (using the calculated instantaneous speed).

My questions are:

Is this approach correct? Or should I just divide the signal into sections corresponding to the cycle of the whole machine and analyse them?
Are there some ready methods to deal with this task, that I could read on?
Can I just cut the signal exactly where the movement of the pinion starts and use it just like that or should I use some windowing function?


texts: signal sampling for 2nd time of sampled signal

texts: A sinusoidal signal of 600hz is sampled with 1khz, If the sampled signal is applied to an ideal low pass filter of 500 hz, what is the output?


texts:
This seems a simple task, yet I am making some stupid mistakes I am not able to find. I create a sinus wave (integer) and then convert it to 16 bit char, little endian. Then I save to file and open it in Audacity (see image below). The first channel in the image is created by Audacity, the second by me. Something is really wrong there. I think it must be in the bit shifting, because if I use 8 Bit, it works. So the extra byte causes trouble but I don't see, how.

#include <bitset>
#include <iostream>
#include <chrono>
#include <math.h>
#include <fstream>
#include <iomanip>

#define PI 3.14159265359

void
test_sine16Bit(){

    int samplesize          = 16;
    unsigned int samplerate = 44100;
    int duration            = 3; // seconds
    int channels            = 1;
    //long samples            = 4096;
    long samples            = samplerate * duration; // per channel
    long total_samples      = samples * channels; // all channels

    // Create test input buffer
    char input[total_samples * samplesize/8];

    // Initialize test buffer with zeros
    for (int i=0;i<total_samples * samplesize/8;++i){
        input[i]=0;
    }

    // Create test input
    int total_index=0;
    char graph[20];
    for (long i=0;i<samples;++i){
        // Channel 1: 440 Hz sinus
        int sinus    = .8*0x8000*sin(2*PI*440*i/samplerate); //  amplitude = 0.8 * max range; max range = 0x8000 = 32768 ( max value for 16 Bit signed int )

        // Little Endian: BUG HERE?
        input[total_index] =  sinus & 0xFF;             // LoByte
        //input[total_index+1] = (sinus & 0xFF00) >> 8; // HiByte
        input[total_index+1] = (sinus >> 8) & 0xFF;     // HiByte
        total_index += 2;

    }

    // Export raw audio data to file
    std::ofstream fout("sine_char16bit.raw");
    if (fout.is_open()){
        long idx = 0;
        for(long i = 0; i < samples; ++i){
            fout << input[idx] << input[idx+1]; idx+=2;
        }
        fout.close();
    }
}

EDIT:

Found the solution: You need to open the file in binary mode and not write via <<, but with put or write.

std::ofstream fout("sine_char16bit.raw", std::ios::out | std::ios::binary);
if (fout.is_open()){
    fout.write(input, total_samples * samplesize/8);
}
fout.close();


texts: Understanding Communication Theory Jargon for Statistical Linear Models

texts:
Coming from a statistical background, I understand the following formula:

\begin{equation}
y = X\beta + \epsilon
\end{equation}

as a linear regression model, with data input-output pair: $(X,y)$, noise term, $\epsilon$, and regression coefficient vector $\beta$.

Now I am confused in how this system is described in communication theory. I see terms such as Rayleigh Fading, MIMO, Multi-access, Multi-antenna, Band-limited etc ... in relation to equations that look exactly like statistical linear regression.

So far I suppose I can work out that for communication theory,

\begin{equation}
y = Hx + \epsilon
\end{equation}

$x$ is some input signal at a transmitter antenna, $\epsilon$ is noise again, $y$ is the signal at the receiver antenna ... but what is $H$? I don't understand how this term should be described, and how it is modelled with respect to all those previous terms I have listed? Reading journal papers I see the simple system as above described with lots of jargon I can't follow, and I'm unaware of the implications this has on the modeling.

If someone could clarify all this jargon and connections for me that would be suuuper helpful!


texts: What is the relationship between the PSD of a continuous signal and the PSD of its periodically sampled one?

texts:
In Oppenheim's "discrete time signal processing", it says,

How is the equation (10.65) derived?
PS:
you may try to derive it with Fourier transform of the original and the sampled signals.
But i think it is not appropriate.
Because when you are considering about the PSD of a signal, the assumption will be the fourier transform of the signal doesn't exist...
the definition of the PSD of a signal actually comes with a truncated signal:

the same is for the discrete version:

I guess there may exist something like... "sampling a random process" theory?
Regards.


texts: Fourier series coefficient of signal when Time period is twice the fundamental period

texts:
My try:

First of all I tried observing the symmetry but I did'nt find any.So I tried to calculate the fourier series coefficient of the signal like this

First I differentiated the signal $x(t)$ so that $y(t)=\frac{dx(t)}{dt}$ Now suppose $d_k$ is the Fourier series coefficient of the differentiated signal and $a_k$ is the fourier series coefficient of Orignal signal $x(t)$
Now $d_k=(jk \omega_0 )c_k$ so $c_k=\frac{d_k}{jk \omega_0}$.Now calculating $d_k$

That is

\begin{align}
d_k &=\frac{1}{3} \int_{-1}^2 \big(\delta(t+1)+\delta(t)-2\delta(t-1)\big) e^{\frac{-jk2\pi t}{T}}dt\\
&=\frac{1}{3}\big(1+e^{\frac{jk2\pi}{3}}-2e^{\frac{-jk2\pi}{3}}\big)
\end{align}

Now $$a_k=\frac{\frac{1}{3}\big(1+e^{\frac{jk2\pi}{3}}-2e^{\frac{-jk2\pi}{3}}\big)}{jk \omega_0}$$
$$a_k=\frac{\big(1+e^{\frac{jk2\pi}{3}}-2e^{\frac{-jk2\pi}{3}} \big)}{j2\pi k}$$

Now checking at

$k=1,a_k\ne 0$
and at $k=2,a_k\ne 0$

so for one odd and one even value of n its not zero so no options matching.What is the mistake?

EDIT:

As stated by @Matt L in the commment to check for period T=6 also, so I did like this

I differentiated the signal from -3 to +3 so the fourier series coefficient I got like this

\begin{align}
d_k &=\frac{1}{6} \int_{-3}^3 \big(\delta(t+3)-2\delta(t+2)+\delta(t+1)+\delta(t)-2\delta(t-1)+\delta(t-2)\big)e^{\frac{-jk2\pi t}{6}}dt\\
&=\frac{1}{6}\big(1+e^{jn\pi}-2e^{\frac{jk2\pi}{3}}+e^{\frac{jn\pi}{3}}-2e^{\frac{-jk\pi}{3}}+e^{\frac{jk2\pi}{3}}\big)
\end{align}

Now at $k=1,d_k=0$ at $k=2,d_k\ne 0$,at $k=3,d_k=0$ but at $k=6$ also $d_k=0$

Now whats the mistake???


texts: Checking the "cleanliness" of an image for Optical Character Recognition

texts:
Is there some mechanism to say whether the image is "clean" for the purpose of image processing?

For example, in the following image, the topmost row is the cleanest and the bottom one is the least clean. As a result, the top-image gives accurate output; the bottom-one gives less accurate.

Examples of noisy and noise-less images

The purpose is to use the image for Optical Character Recognition.
And, I wanted to check beforehand, if the image is sufficiently clean or not.


texts: Is it possible to recover a waveform from spectrograms of magnitude *and* phase?

texts:
It is my rough understanding that a typical spectrogram discards phase information by squaring only the real-valued part (magnitude) of the STFT (I know it's not quite this simple, but feel free to correct me if I'm wildly wrong here).

While many software libraries (e.g., scipy) will, by default, produce a spectrogram of magnitude, they often also make it possible to instead compose a separate spectrogram for phase. My question is quite simple: if I had both of these matrices, could I reconstruct the original signal?

I suspect the answer to this question is 'not quite'. If so, what part of the process is problematic and are there any changes that could be made to make this goal achievable? I am currently using STFT.real and STFT.imag (where STFT is a 2D matrix for the short time Fourier transform of some signal) and plotting of them separately. While I can reconstruct the original signal using this approach (which I need to be able to do), these plots are not as easy to read as actual spectrograms.

(For instance, I know that squaring the magnitude is clearly a problem for reconstruction as it destroys negative numbers).


texts: Analog vs. digital communications in PCM and noise

texts: How do you determine if you are using analog or digital communications? For example when dealing with PCM and noise how can one notice any difference in the signals when using the two types of communications? Thanks


texts:
I have read many of the articles on the website regarding a polyphase decimating filter structure.  I'm sure we are all familiar with the typical image shown in Fred Harris Multirate book and many other papers written on this topic below.

Say we have some transfer function H(z).. well the polyphase implementation works by decimating the filter into M parallel filters where fs/M is the desired output rate.  The noble identities tell us we can pull the decimation through the filter stages... and hence the delays at the input are replaced by the commutator switch which effectively sends samples in at rate fs/M.  Let N be the total number of taps in H(z)... the convenience of this is that we only need to filter with N/M taps per cycle which greatly reduces resources... ok I've done my homework.

Now my question.  I have seen many implementations of this structure in code but I can't seem to find any implementations in simulink...specifically with what exactly is the commutator switch?

Here is an implementation that I have created in simulink... although it clearly doesn't work properly...Say I have some filter that was originally 120 taps and I now want an output rate that is 1/4 of the input rate... I decimate my filter into 4 sub filters and arrange into a parallel structure as shown below:

The switch on the input has a counter attached to it because I'm thinking that I can cycle between inputs at each clock edge to simulate the commutator...but then I need to demux the data into a branch that I want.. however this isn't working because the demux isn't controllable...

What am I doing wrong here?  What is the block that simulates the commutator like shown in many of the block diagrams in the multirate papers we all know?  I realize that it's a memoryless switch but which simulink block ?


texts:
The design of control systems, particularly for SISO systems, is made convenient by tools such as Bode plots and the corresponding Nyquist stability criterion. These concepts allow engineers to design sufficient margins for the given system in a model based paradigm, when it's transfer function is known. Specifically in the Bode plot, the frequencies correspond to sinusoidal input signals and the plot provides information about the gain and phase margin available for stability. The idea is that by knowing the worst case of frequencies one can design the controller to reject disturbances.
However, in the case of stochastic noise, say the commonly used model of Gaussian white noise, it is not clear to me what the corresponding frequency representing this type of noise would be on the Bode plot. What is the frequency of the Gaussian white noise model in control systems? The first variation (time derivative) of this noise is infinite indicating very large frequency. Additionally, does this fact point to the evolution of stochastic control as a separate sub-discipline within robust control (control methods to reject adverse effects of disturbances and noise)? Finally, is there a frequency-domain representation of Gaussian white noise?


texts: Creating Bode Plot from Experimental Data

texts:
I have a blackbox system in which I can input a function and obtain the output signal (in MATLAB). I'm attempting to reconstruct a Bode diagram and have had success with the Magnitude plot, however, I'm having a lot of difficulty with the phase plot. What I've tried to do is collect input and output data for the input $\sin(ft)$ where $f$ is the frequency that I am changing. For each frequency I am determining the phase shift by finding the time difference (I'll denote by $td$) between the input and output signal on a plot and calculating the phase shift using
$$td \div\frac{2\pi}{f}$$
After collecting the data over a range of frequencies I plot this (phase shift on y-axis, frequency on a log scale on the x axis). The plot does not at all look like a phase plot, I don't really see much of a pattern, nor does the scale seem correct for the phase shifts. Below are some examples of data I have collected for different frequencies:
Frequency: $f=0.05$: https://ibb.co/swCGLw8
Frequency: $f=10$: https://ibb.co/jL8JhTv
Frequency: $f=1000$: https://ibb.co/rkZ8WLh
Note that a filter has been applied to try and reduce the noise
Any help would be great!
EDIT: Image of the phase plot: https://ibb.co/K2LcHqJ


texts: Confusion about result of FIR All-pass filter design

texts:
I try to design a FIR All-pass filter with random phase in the frequency domain. I am a bit confused by my result and am not sure if the reason is a programming error or a misconception about convolution and impulse response of myself. I use Python and Numpy:

import numpy as np

mag = np.ones(129) # magnitude is unitiy
phase = np.zeros(129)
random_numbers = (np.random.rand(128)-0.5)*2.0
phase[1:] = random_numbers * np.pi # phase is random except for dc offset
spec = mag*np.exp(phase * 1j) # get spectrum from magnitude and phase
coefficients = np.fft.irfft(spec) # do inverse fft to get coefficents, result is of length 256

If I do a fft on my coefficients the magnitude looks nearly right. My first question: Why is the highest frequency suppressed?

And secondly, if I use my coefficients for convolving with an impulse, I thought that should reproduce the coefficients? But it does not. Why are there these "random variations" around 1.0?

impulse = np.zeros(256, dtype=np.float64) # create impulse
impulse[256/2] = 1.0
result = np.convolve(impulse, coefficients, "same") # convolve it

Any explanations or hints are much appreciated!


texts:
I tried to read up and get a filter for a 3-axis accelerometer but this is more confusing than I thought.
I have an IMU with a 3-axis accelerometer on the end of a robot link which is rotated  around on one axis.
When at rest with just gravity acting on it I get these readings:

Now the band of the data is kinda broad and very noisy . I thought that at first I could maybe look at the frequency spectrum and cut out unwanted frequencies but honestly the FFT doesn't tell me anything I could remove. What does initial spike in the spectrum mean and and how would I put a high or low pass filter over the data if it looks like such a mess? Do I need to smooth it in some way?
Preferable I would want a nice arc of the data which follows the movement of the end of the arm but I don't know how to achieve that. Could it be that the arm movement by the motor introduces too much noise?
Thanks for any help!


texts: Transform minimum phase FIR into linear phase FIR

texts:
I've seen examples of transforming a linear phase FIR into a minimum phase FIR, but is there a simple process to transform a minimum phase FIR into a linear phase FIR?

I would like to end up with a single FIR of linear phase. If the magnitude ends up being squared during the process (or the number of taps double), then that would be OK.


texts: Why Does the Rudin Osher Fatemi (ROF) Method Use Variational Methods for Image Denoising When Denoising Problems Are Not Boundary Value Problems?

texts:
I was reading the recent book Variational Methods in Image Processing by Vese and Guyader which is quite interesting. In the book the authors discuss different types of image processing problems, from denoising and deblurring problems to image segmentation and beyond. The common theme is that they use Variational Methods, read as Calculus of Variations methods to take the noisy image and remove the noise. My question is, why do these scholars use Variational Methods (Euler-Lagrange Equation) to solve a problem that is not a Boundary Value Problem. Let me lay some foundation to the question.
Now these image problems are naturally optimization problems where I have some original image $f$ and I want to find a new image $u$ that best replicates $f$ but is subject to some constraints. The Rudin-Osher-Fatemi method sets up the problem this way:
$$
F(u)=\int_{\Omega}|D u|+\lambda \int_{\Omega}(K u-f)^{2} d x
$$
Where the second term is the standard difference between the original image $f$ and the adjusted image $Ku$, and the first term is the Total Variation norm in the image, meaning the sum of the absolute value of the first difference between adjacent pixels.
The solution to this problem sets up the Euler-Lagrange equation for the above energy function, and then minimizes that energy using gradient descent or such algorithms.
BUT, in my mind Calculus of Variations was always taught as a way to solve a boundary value problem. I am given two points, and I want to find the least energy path between those two points. So I have a BVP and I can either solve the optimization problem directly by discretizing the space between the two points and using shooting methods + opimizer to find the minimal energy path. Or I would use the Euler-Lagrange equation to find the least energy path as the numerical solution ODE or PDE--and thus avoiding minimization.
In looking at image processing problems, I did not have a sense of a boundary value problem? The initial image is fine but the final image is the result of optimization and not known before hand. Hence, can anyone explain why is the Variational approach used, or Euler-Lagrange approach?
Are they just using the E-L to find a minimum without thinking of a boundary value, or is there some other intuition going on that I am just not catching?


texts: How to sketch a frequency spectrum for an AM signal?

texts:
So I have an AM signal, given by:
$$v(t)=50\cos(\pi 10^{6} t )+ 20\sin(\pi 10^3t)\cos(\pi 10^6 t)$$
I was asked to sketch the spectrum of the signal obtained, so I found the fourier transform of it, and that's what I got,

$$V(f)=25 \left[ \delta(f - 500k)+ \delta(f+500k)\right]-5i\left[\delta(f-999.5k)-\delta(f+999.5k) \right]-5i\left[\delta(f-500.5k)-\delta(f+500.5k) \right]$$

And so far I only knew how to sketch the first term, as shown below,

But I don't know how to sketch the terms with the $-5i$.


texts: Telling apart a right traveling non-periodic signal from a left-traveling non-periodic signal, with a background periodic signal

texts:
I have a time-signal sampled over a one-dimensional spatial domain $x$. I know it is made of 3 components:

a right traveling non-periodic wave (characteristic time width $T_{R,NP}$, mostly 0 outside this width). This is prescribed at the left end of the domain.
a left-traveling non-periodic wave (characteristic time width $T_{L,NP}$, mostly 0 outside this width). This is prescribed at the right end of the domain
a periodic signal (right traveling in my case, of period $T_P$)

While traveling, all these signals slightly change shape due to non-linearities. Again changes are not drastic and fairly small. The periodic signal starts as a single harmonic but then more harmonics show up, which have small energy and the peak period is still $T_P$. Here is an example of the signal at a specific location (purple), together with my best guesses (I drew them by hand) of the two non-periodic signals. Two signals were prescribed on the left and traveled right, and one was prescribed on the right and travel left.

Which is the best way to tell apart the 3 components at different locations of my domain? In particular:

Is there a method that can take advantage of the known information about the direction of travel?
Sometimes (but not always) $T_{L,NP} > T_{R,NP} > T_P$. Would this help? I would prefer a general algorithm that does not rely on this information, but if a general algorithm does not exist, solving for these special cases would still be something.


texts: Should scipy.signal.coherence be 1 for single input and output signals?

texts:
I am trying to calculate the coherence between input and output signals. I thought I could work with a single input and a single output time series and calculate the coherence $\gamma^2$ between them. I did so with the scipy.signal.coherence function (See code and results below).
Now, a colleague challanged my procedure, saying that the coherence needs to be calculated by averaging several input and output signals, since otherwise the coherence will simply be $\gamma^2=1$. I wanted to prove otherwise, after all I used a scipy function and I regard those as trustworthy.
So I started out with the general definition of the coherence as [1][2] [3]:
$$
    \gamma^2 = \frac{|\mathbb{E}[X(-j\omega)Y(j\omega)]|^2}{\mathbb{E}[|X(j\omega)|^2]~~\mathbb{E}[|Y(j\omega)|^2]}, \quad \hat{\gamma}^2 = \frac{|\overline{S}_{XY}|^2}{\overline{S}_{XX}\overline{S}_{YY}},
$$
where $\overline{\cdot}$ indicates a mean value, $|\cdot|$ the magnitude of a complex number, ${S}_{XX}$ and ${S}_{YY}$ are the energy spectral density for $X$ and $Y$ respectively and ${S}_{XY}$ is the cross spectral denisty between $X$ and $Y$. And $X(j\omega)$ being the Fourier Transform of a time series $X(t)$ and $X(-j\omega)$ being its complex conjugate. Let's define those as:

$$
X(j\omega) = x_1+ix_2, \quad X(-j\omega) = x_1 -ix_2, \quad Y(j\omega) = y_1+iy_2.
$$

Now, assuming I do not have several signals, but only a single input signal $X$ and single output signal $Y$, I drop the mean $\overline{\cdot}$ and end up with the following:
$$
    \gamma^2 = \frac{|{S}_{XY}|^2}{{S}_{XX}{S}_{YY}} = \frac{|X(-j\omega)Y(j\omega)|^2}{|X(j\omega)|^2|Y(j\omega)|^2} = \frac{\left(|X(-j\omega)||Y(j\omega)|\right)^2}{|X(j\omega)|^2|Y(j\omega)|^2}= \frac{\left(\sqrt{x^2_1+x^2_2}\sqrt{y^2_1+y^2_2}\right)^2}{\left(x^2_1+x^2_2\right)\left(y^2_1+y^2_2\right)}=1.
$$

In case I did not make a mistake, my scipy.signal.coherence function should give me a coherence = 1. But as shown below, it does not. Instead the coherence is close to zero everywhere, except at the signals' natural frequencies. Please note that I get $\gamma^2=1$ only if I chose $X=Y$.

I must be missing something here but I just don't see it. Any ideas?

Many thanks in advance!
Max

Result:

Code:

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
from scipy.fftpack import fft

fs = 5e3 # sampling frequency
N = 10000 # number of observations
T = N/fs # max time T
amp = 20 # sine wave amplitude input
amp2 = 2*amp # sine wave amplitude output
freq = 500.0 # sine wave frequency input
freq2 = 500 # sine wave frequency output
time = np.arange(N) / fs # time vector t_0 up to T
# create data
x = amp*np.sin(2*np.pi*freq*time) # input
y = amp2*np.sin(2*np.pi*freq2*time)+np.random.normal(scale = 10, size = len(x)) + 100 # output

fig, axes = plt.subplots(1,3, figsize=(15,5))

# plot time signals
axes[0].plot(x[:400], label = 'Input', alpha = 0.7)
axes[0].plot(y[:400], label = 'Output', alpha = 0.7)
axes[0].legend()
axes[0].grid(linestyle='--')
axes[0].set_xlabel('Time')

# plot fft input
xfft = fft(x)
N = len(x) # number of sample points
T = 1 / fs # sample spacing
xfreq = np.linspace(0.0, 1.0/(2.0*T), N//2)
axes[1].plot(xfreq, 2.0/N * np.abs(xfft[0:N//2]), label = 'Input', alpha = 0.7)

# plot fft output
yfft = fft(y)
N = len(y) # number of sample points
T = 1 / fs # sample spacing
yfreq = np.linspace(0.0, 1.0/(2.0*T), N//2)
axes[1].plot(xfreq, 2.0/N * np.abs(yfft[0:N//2]), label ='Output', alpha = 0.7)
axes[1].legend()
axes[1].set_xlim(1,1000)
axes[1].grid(linestyle='--')
axes[1].set_xlabel('Frequency [Hz]')

# plot coherence:
f, Cxy = signal.coherence(x, y, fs)
axes[2].semilogy(f, Cxy, color='black')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Coherence')
axes[2].grid(linestyle='--')

plt.show()


texts:
I am trying to compute the fundamental phasor using sliding window DFT. I have employed a Blackman window in conjunction i.e
$$
\sum_{k=0}^{L_{DFT}-1}x(k) w(k) e^{-j2\pi k/N}
$$
where $x(k)$ is the time-domain signal sampled at a rate which is an integral multiple of the fundamental frequency i.e.$F_s$=$(f_0 \cdot N)$, where $N$ is the number of samples per cycle of the fundamental, and $w(k)$ represents the Blackman window.
I am referring to a paper which is giving me an equation that should compensate for the magnitude attenuation due to this window whenever the time-domain signal frequency deviates from the nominal. The equation is
$$H_{2}(e^{j \omega_2})=0.42 H_{1}e^{j\omega_2}
+0.25\left(H_{1}e^{j\left(\omega_{2}-\frac{2\pi}{L_{DFT}-1}\right)}+H_{1}e^{j\left(\omega_{2}+\frac{2\pi}{L_{DFT}-1}\right)}\right)
+0.04\left(H_{1}e^{j\left(\omega_{2}-\frac{4\pi}{L_{DFT}-1}\right)}+H_{1}e^{j\left(\omega_{2}+\frac{4\pi}{L_{DFT}-1}\right)}\right).$$
The author also mentions that
$$
\left\lvert H_{1} e^{(j \omega_{1})}\right\rvert=\left\lvert\frac{\sin(L_{maf}\omega_{1}/2)}{L_{maf}\sin(\omega_{1}/2)}\right\rvert$$
where $L_{maf}$ is the length of a simple moving average filter used in the pre-DFT computing process and
$$
\omega_1=\frac{2\pi f_F}{F_s L_{maf}}\quad\text{also}\quad \omega_{2}=\frac{2\pi f_{F}}{F_s}$$
with $f_F$ being the computed time-domain signal frequency.
Using values such as $L_{maf}=31$, $L_{DFT}=95$, $f_0=50$, $F_s=50\cdot 24$.
The list of values obtained for $H_2$ for a list of values of $f_{F}$ are given below:

50 Hz=0.0639,
51 Hz=0.0618,
52 Hz=0.0594,
53 Hz=0.0569,
54 Hz=0.0543,
55 Hz=0.0516.

Questions:

Why is $L_{maf}$ involved in $H_2$ computation?
Why are the values of magnitude gain of the blackman window not accurately compensating for the magnitude attenuation caused by the same, i.e. for 50 Hz?
The gain should have been unity exactly since the DFT is centered around it, and then close to unity as the frequency deviates from the nominal.


texts: How do I find the ROC of a system if it has no poles

texts:
The output of a system of discrete time $y[n]$ is corellated with the input $x[n]$ through the equation $y[n]$.

$$y[n] = \frac 13\big(x[n-1]+x[n]+x[n+1]\big)$$

It then asks me to find the system function, the ROC, the zeroes and the poles.

I found out that $$Y(z) = \frac 13\left(z^{-1}X(z)+X(z)+zX(z)\right)$$

I found out that the system function is : $$H(z) = \frac {Y(z)}{X(z)} = \frac 1{3z}(z^{2}+z+1)$$

The pole here is 0. How do I tell if this sequence is left/right/two sided to figure out the ROC?

This is an unsolved exercise given by my professor.


texts: Spectrum of WBFM signal for triangular message signal spectrum

texts:
We Know:
Frequency Modulated signal is given by:
$$s_{FM} (t)=A_{c} \sin \{ w_c t + 2 \pi k_f \int m(t) \, dt \}$$
Now,
if $ \quad m(t)= A_m \cos(w_m t)$ ,
i.e, its spectrum $M(f)$:

then,
$$s_{FM} (t)=A_{c} \sin \{ w_c t + 2 \pi k_f \int A_m \cos(w_m t ) \, dt \}$$
$$= A_{c} \sin \{ w_c t + 2 \pi k_f \frac{ A_m \sin(w_m t ) }{w_m} \}$$
$$= A_{c} \sin \{ w_c t + \frac{ k_f  A_m  }{f_m}  \sin(w_m t ) \}$$
$$=A_{c} \sin \{ w_c t + m_f  \sin(w_m t ) \}$$
$$= J_0(m_f) A_c \sin(w_c t) + \sum_{n= 1}^{\infty} [ \{  J_n(m_f) A_c  \sin(w_c t + n w_m t) \} + \{ J_{-n}(m_f) A_c\sin(w_c t - n w_m t)\}]$$
$\implies$ the corresponding FM (Frequency Modulation) signal spectrum [$S(f)$] will be:

Now, if the spectrum of message signal is triangular,
i.e,

then, what will be its corresponding WBFM (Wide Band Frequency Modulation) signal spectrum?
i.e, $m_f > 1$

for e.g:
For $m_f \leq 1$ , the corresponding NBFM signal spectrum is:


texts: Pole zero plot, normalizing frequency response plot?

texts:
I'm asked to plot the frequency response (amplitude) given a specific pole-zero diagram.

$$ H(z) = H_0 \frac{\prod\limits_{m=1}^{M} (z - q_m)}{\prod\limits_{m=1}^{M} (z - p_m)}$$

$$ H(e^{i\omega}) = H_0 \frac{\prod\limits_{m=1}^{M}(e^{i\omega} - q_m)}{\prod\limits_{m=1}^{M}(e^{i\omega} - p_m)}$$

If I understood it correctly, the amplitude at frequency $\omega$ is (the magnitude of the distance from $e^{i\omega}$ to all the zeroes) divided by (the magnitude of the distance from $e^{i\omega}$ to all the poles), i.e:

$$ \Big|H(e^{i\omega})\Big| = \Big|H_0\Big| \frac{\prod\limits_{m=1}^{M}|e^{i\omega} - q_m|}{\prod\limits_{m=1}^{M}|e^{i\omega} - p_m|}$$
where $q_m$ are the zeroes and $p_m$ are the poles and $H_0$ is the constant gain factor.

The problem is that in the graph where I need to draw the frequency response, the frequency and amplitude range $0\to1$ like so:

After I get a value by calculating the poles and zeroes I almost always get a value above $1$. What do I need to do with the value to fit it into the graph? How do I normalize(?) the value?


texts: How can I start learning about pitch detection?

texts:
readers.

I'm working on a project where my goal is to implement a pitch-detection algorithm (PDA) in a piece of open-source software.

I have very little knowledge of pitch detection at the moment and started my research into it by reading through part of a paper on PDAs.

I wrote a similar post to this one on Music: Practice & Theory Stack Exchange (SE) and the advice I got was to start learning about digital signal processing.

My current plan is to study a set of lecture notes I have on digital signal processing.

I was also directed from Music SE to this SE site. Would any of you suggest a different or potentially better approach to learning about digital signal processing and PDAs (keeping in mind my ultimate aim of implementing my own PDA - I don't want to simply copy or plug in a library of someone else's, but study the pros and cons of other PDAs and then implement my own, with some element of novelty)?


texts: Computing the center of mass of a signal using the Fourier transform

texts:
I am having difficulty with a homework problem which asks:

I have some ideas in mind but I have no clue as to whether they are correct or not. Below is my attempt:


texts:
I am after a simple, yet effective, Image Edge Preserving Filter.
I need something which is faster than the Bilateral Filter with similar quality (The Guided Filter isn't good enough).
Are there such algorithms?


texts: Difference between a wavelet transform and a wavelet decomposition

texts:
I'm confused about the difference between a wavelet transform and a wavelet decomposition is. For example

load woman
[cA1,cH1,cV1,cD1] = dwt2(X,'db1');
[c,s] = wavedec2(X,2,'db1');

What's the difference between these two matlab commands, and when would you want to do one over the other?


texts: Image Subtraction using Pre-bayer data (.pgm)

texts:
Not sure if this is the best stack exchange but I'll give it a shot.

I am performing background subtraction to collect Electroluminescent emission (i.e image w/ emission - image w/o emission = emission). Anyway, I am using a DSLR and collecting these images in .cr2 (canon raw format). Since the emission is weak, we wan't to make sure that we have the most data available to us and that it is as unprocessed as possible so we are converting the RAW images to .pgm (pre-bayer image data).

There is an issue, in some of the background that does not emit EL, for example a cabinet with drawers, when we perform background subtraction, we can detect all of the edges. Additionally, even if I try to do subtraction on 2 background images (both have no emission), the edges can still be clearly detected where we would expect simply black areas in the image.

Would anyone have insight on whether this is due to using .pgm, etc. Any suggestions at all?

Feel free to ask more questions about the project etc.
Thanks

Below is what I get subtracting 2 images which should be essentially the same. I capture the raw image .cr2, convert it to .pgm using dcraw then subtract them using python numpy arrays and output a viewable png (what you see below.). Although, we would expect a black image(i.e subtracting two same/similar images), here I can detect all of the edges which makes no sense to me. Any suggestions?


texts: Why is it used the Power Density Spectrum (PSD) over an anylisis with the Fast Fourier Transform (FFT)?

texts: I'm currently working with physiological signals (PPG and GSR) for emotion recognition but, from my research, I've found out that almost everyone in that area use a PSD analysis over a FFT analysis. I've been reading about them and found out that PSD helps with giving a clearer view of the spectrum despite the amount of data that you have, as described in this blog https://blog.endaq.com/why-the-power-spectral-density-psd-is-the-gold-standard-of-vibration-analysis, also because it is supposed to be used with random signals and therefore GSR and PPG are signals that have some random nature within them. Despite this, I still can't gasp the intuition why a PSD is used over a FFT analysis.


texts:
I have been experimenting some weeks to find a way to match/find similar Songs in a Library containing different Genres of Music.

My first try was to Detect Features like Tempo or how much Bass there is in Songs to form groups, but I didnt get far with this approach (Volume Changes Based Beat Detection) since on about 20% of the Songs Beat dont have to be counted always, sometimes 1/2 or 1/3 of them and I couldnt implement that.

After some weeks of failed trying I got a new idea which is described later in this Post. Simply put it works by taking Spectrum Samples of Files, making something like an "Average Spectrum" of Files to compare them.
The Idea behind was that for example Hardstyle has a lot more Bass than average Rock Music, I also verified this by looking at some Spectrums in Audacity.

File 1: Take full File FFT Spectrums (2048 Sample Size atm,Amplitudes Log. scaled)
Sum all Spectrum Arrays, take Averages of each Bin
Do the same some other Files, store all Results
Make List of FFT Values Differencies between File 1 and other Files
Make Average of Differencies between File 1 and File X
Sort Ascending by these Averages
Songs with the Lowest "Difference Value" are considered to be Similar.

Can some of you who have good knowledge tell me if this would be the right/good way to implement my Idea?


texts: Filter wrap-around to avoid dealing with settling time? Using a Butterworth

texts:
I've been using a Butterworth high pass filter to correct the low frequency distortion of an acceleration waveform integrated to a velocity waveform.  I throw out a small number of the first data points of the velocity waveform due to the filter settling time.

Is it acceptable or even possible to wrap the filter back around to the beginning and process those 'settling' points as well so I don't have to throw out any data?

Maybe you can reverse the direction of the filter when you hit your known settling point and process back to the beginning then reverse again and go all the way to the end?

I'm not sure if this is an obvious or appropriate question, I have no background in DSP.


texts: How can we measure the quality of the FIR filter?

texts:
The process of designing filters includes too many parameters like bandwidth, transition, degree of attenuation in side-bands, phase delay, causality, the magnitude of ripples, and so on. First of all, we want to approach the ideal filter structure so we can design many different filters to obtain the desired ideal filter.
How can we compare all of these filters in terms of quality? Is there any way to measure the performance of filters? I am looking for a mathematical solution to this problem.


texts: Count the vehicles flow that pass over one bridge

texts:
I live in an apartment near a bridge. From my window I can see the cars traveling over the bridge as you can see below:

Is there a relatively simple manner to count each car pass over the bridge? Using an android cell phone pointed to the bridge behind a lunette and using a API for image processing like Open CV?


texts: How do I identify a notch filter given the transfer function?

texts:
I am studying for an exam, and a problem that my professor gave produces a transfer function that looks like the following:

$$H(z)=\frac{z^2-j}{z^2-\frac{1}{4}j} = \frac{\left(z-e^{j\frac{\pi}{4}}\right)\left(z+e^{j\frac{\pi}{4}}\right)}{\left(z-\frac{1}{2}e^{j\frac{\pi}{4}}\right)\left(z+\frac{1}{2}e^{j\frac{\pi}{4}}\right)}
$$

How am I able to determine:

That this is a notch filter,
At which frequencies the notches occur at?

His solution says that the notches occur at $\frac{\pi}{4}$ and $-\frac{3\pi }{4}$.


texts:
How can I use interpolation with an FIR filter? I am more familiar with interpolation such as nearest distance interpolation, linear interpolation and so on.

Suppose a signal is given as the following series:

    time=    1 2 3 4 5 6 7 8
    signal=[ 1 1 2 2 3 2 7 6]

Then, if I ask what's the signal value at time 2, we can see signal[2] = 1; then how about the signal value at time 2.1? With a linear interpolation method, I can tell signal[2.1] = 0.9*1+0.1*2 = 1.1.

In this case, how could I know the corresponding FIR filter? How could I use the filter to obtain the signal value?


texts: Image processing (number recognition) woes

texts:
So I have the following image:

Using GIMP, I was able to get to here:

Select all white color (with a touch of fuzz)
Invert selection
Cut selection
Invert all color (goes black)

But now I'm stuck...

I want to use something like ocrad to extract the numbers.

ocrad goes crazy when I input the black&white image above.

I'd like to see (the grammar is the numbers 1 to 45 inclusive):

var a = [30,31,32,35,37,40,44];
var b = [6,7,11,15,18,21,22];
var c = [5,11,15,18,23,37,28];

It's worth pointing out that the image

is always the same dimensions
is always the same orientation
therefore, I think I can split the image into three rows relatively robustly

Like so:

(in the above case, I also did an extra step - remove black contiguous regions that have large area. Though I don't know how to do this programatically - I did this by eye in GIMP using the "magic wand" tool).

The above image does OK in ocrad, but it gets into trouble with the ring surrounding the "44".

I'm trying to avoid pattern matching.

Does anyone have any tips? I'm thinking of doing all this in Python/OpenCV/similar.


texts: Identifiability of a state space model (Dynamic Linear Model)

texts:
Take a general linear Gaussian state space model (SSM)(aka Dynamic Linear Model DLM):
\begin{align}
X_{t+1}&=FX_t + V_t\\
Y&=HX_t+W_t\\
V_t &\sim \mathcal N(0,Q)\\
W_t &\sim \mathcal N(0,R)
\end{align}

I am interested in the unidentifiability issues related to these models:

Hamilton (1994) states that:

  in the absence of restrictions on $F$, $H$, $Q$ and $R$, the parameters of the state-space representation are unidentified - more than one set of
  values for the parameters can give rise to the identical value of the
  likelihood function, and the data give us no guide for choosing among
  these

Now I realise that this representation is not unique as multiplication by any orthonormal matrix $M$ produces a new representation:

\begin{align}
MX_{t+1}&=MFM^{-1}X_t + MV_t\\
Y&=HM^{-1}MX_t+W_t
\end{align}

This type of unidentifiability where the observed values can be produced by various orthonormal transformations of the state variables is inherent to state space models.

However, I also came across another type of unidentificability that seems to be related to the estimation method. In this case "Kalman Filtering". See the simple example starting on page 8 here.

In this case there is a linear transformation of the observation equation and an offsetting one made to the variance of the state equation

Do both the transformations above give rise to the same kind of identifiability issues described by Hamilton (I believe they do but want to check) ?
Are there other ways in which identifiability issues can manifest themselves in Linear Gaussian SSMs ?
Is the fix always the same find constraints or analogously (Bayesian priors) that will ensure that the final parameters are correct ?

Lastly, this link in Matlab suggests that it is possible to build an "identifiable SSM":

Unfortunately the link doesn't explain the theory:

It is possible to translate any Linear Gaussian SSM into an "identifiable form" ? Can someone please supply a link reference explaining how this works. It would seem at first blush that no matter what representation is used initially it would still be subject to the problems shown above ?

Thanks

Baz


texts:
What are disadvantages of root raised cosine pulse shaping filter in digital communications and why does it need to be improved?

Links:

Square Root Raised Cosine Fractionally Delaying Nyquist Filter
Square Root Raised Cosine Filter
Generalized raised-cosine filters


texts: Why are there vertical lines in the power spectrum?

texts:
I took an image, did an fft2 to it in Matlab and then squared the abs of the fourier transform to get the power spectrum.

Something I notice about the power spectrum however is that there are vertical lines in it like this:

This is the original image which the power spectrum was calculated from, for reference:

What are these lines and why do I see such phenomenon from the transformed image? Do the lines imply anything about the original image and the fourier transform?


texts: Non-uniform resampling - spacing between output samples given speed factors

texts:
I'm working on an implementation of time-varying resampling in python. I've got the windowed sinc interpolation itself working nicely, but can't figure out where exactly to sample my input data for correct results.

Given a curve of speed factors speeds (where 1.0 means no speed change, 2.0 slows down and .5 speeds up the result), I interpolate linearly between speeds[i], speeds[i+1] and then take the cumulative sum of the reciprocal of the result to get my positions.

block_speeds = np.interp(np.linspace(0, 1, num_output_samples), (0, 1),(speeds[i], speeds[i+1])  )
positions = np.cumsum(1/block_speeds)

My ground truth assumption using a different resampling method is repeating every input sample (N * speed) times, and then downsampling everything by N.

The approach outlined above approximates my ground truth for short segments (say 0.01 seconds), but strays from it too much for longer segments (say 1 second).

Any idea how to determine the sample positions for the output more accurately and analytically?


texts: How to quantify area of interest in an image?

texts:
How to find the percentage of a certain colour in a image accurately? For example  I am actually trying to find the percentage of barnacles in this image over the red part of the hull.
http://thumb1.shutterstock.com/display_pic_with_logo/1498544/150482510/stock-photo-barnacles-barnacles-perched-on-the-red-buoy-150482510.jpg


texts: why the number of filter coefficients in FIR filter has to be an odd number?

texts:
I calculated the order of the FIR filter to be 31(N). so, the number of coefficients has to be 32 (N+1). so, I have to increase it to 33 to make it an odd number.

Why the number of filter coefficients is required to be an odd number?


texts: Compensate for residual resampling error with a quadratic phase vs frequency- possible?

texts:
With regards to the answer in this post Is it possible to do deconvolution with two data sets that have different sampling rates?

I suggest the possibility of compensating for a residual resampling error with a quadratic phase versus frequency function. I haven't seen this done and have not confirmed its effectiveness or the limits of its use. My question is if anyone has done this or seen it used (as if possible I would imagine it to be commonly used), or if there is clarity on why such an approach would be flawed.

The considered approach would be to remove residual frequency sampling error with use of a all-pass function that has a quadratic phase versus frequency response. (A linear phase versus frequency is a fixed time delay, so a quadratic phase versus frequency would have a time delay that linearly increases versus frequency, in such a direction to compensate for the phase error that linearly increases or decreases versus frequency depending on the sign of the error).

I realize that it may not provide a perfect compensation since the result would be a linear delay versus frequency, but am imagining utility to the same extent that we can approximate a phase shift over narrow frequency ranges with a fixed time delay (it is the desired phase shift at one frequency and then deviates linearly from there). In the case of applying this to a fixed length FFT, note that the effect of a shift in the sampling rate will cause each frequency bin to shift by a certain phase (quadratically dependent on the frequency), so the compensation at each frequency bin should be exact.

I feel I may be missing something obvious. Could this approach reduce residual frequency sampling error?


texts:
I am writing a C++ program to use to analyse oscillations and create a power spectrum.
I am using FFTW3.

To test the system, I am using a perfect sine function.

But no matter what parameters I try, I run into the same problem. I find the correct frequency, but also many 'overtones'.

First I assumed the problem was my way of chopping up the signal or how I calculate the running mean. But the problem is there when I just have all data in a single big window, bypassing my suspect code, still creating the pattern in the pictures.

I then assumed it was an aliasing problem, but even 1024 datapoints in a period, and they don't go away. They just lose shape and appear randomly. I assumed it was too little periods in my window. But more periods per window and they just become less frequent and stronger.
Then I assumed it may be spectral leakage, so I assumed I may need to use a windowing function like a Hamming, but that creates an odd slant. Is that the trade-off?

Or do I just set everything that is small enough to zero? I was just expecting a single strong point peak and the remainder to be so near zero, the log is -6 or lower everywhere else uniformly.

If I can't figure out how to get the right parameters for the best signal-to-noise ratio, I am not sure how I will achieve that when I start to work with the actual data.

I assumed it would be simple to pick a window that contains 32 periods with 32 discrete points in each period, use Welch's method and be done with it.

I am not an engineer who took many many classes on this kind of stuff. So I just don't have the proper training. I am working on a project of 4 weeks. There is so much info out there, it is confusing, but I don't know where to start as I don't have the time to learn stuff properly.

It can't be caused by the limited numbers of decimals of the input sine values, right?


texts: Difference(relation?) between filter banks and wavelet decomposition

texts:
I'm implementing multirate system and met concepts of filter banks and polyphase decomposition of filters.

While trying to understand them, I've found that they look very similar to wavelet decomposition, and now I'd like to understand their relation and differences.

For me, the only thing is kernel types of filters (usually we use sinc kernels in filters, and other, like Haar, in wavelets).


texts:
I have written a program that uses sinc interpolation to resample some data. The general algorithm is a that I compute the previous N values and the next N values to get a new sample at a non-integer value.

For example if I have samples:

(0, 8), (1, 5), (2, 3), (3,8), (4,-2), (5, 11)

And I want a sample value at 2.3, I offset the 2.3 incrementally by 1 and take the sinc coefficients and weight them by the sinc function value e.g.:

sinc(-0.7) * 8, sinc(0.3) * 5 + sinc(1.3) * 3 + sinc(2.3) * 8 + sinc(3.3) * -2 + sinc(4.3) * 11

This algorithm fits my understanding of the sinc interpolation concept, and I'm pretty sure there are no bugs.

Unfortunately I'm hear some tin-y artifacts in the sound.

Given the above algorithm is working properly and is a valid way to resample data (a big premise), is there any reason sinc interpolation would create these artifacts? When I use the same method to down sample by a large enough factor the artifacts to go away.


texts:
I am trying to perform a of convolution 2 3-dimensional vectors using the FFTW3 library in C.

My question relates to the order in which the frequency domain appears and if it will make a difference to the final result.

In the FFTW documentation it dictates that in 1D complex FFT's the negative frequencies are in the second half of the array. Does this follow through with a 3D complex FFT as the documentation doesn't say? If this is the case how would it be stored? i.e in 1D it would be something like:

|+|+|+|+|-|-|-|-|

so what would it be in 3D (sorry, terrible as visualising these things).

Finally when convolving 2 3D vectors using the convolution theorem do I need to do anything because of the FFTW layout, or can I just do the pointwise product and then IFFT back?


texts:
I have to do a project for college. The project consists in writing a paper about some application of digital signal processing and then implement the solution in MATLAB. It doesn't have to be the original.

We're supposed to check papers to come with ideas but I thought that perhaps someone here could give some ideas.


texts:
I have designed a simple 3-state $[x\ \dot{x}\ \ddot{x}]$ Kalman filter which is updated by measurements for $x$ and $\dot{x}$. The filter is tracking a peak on a surface in the $x/\dot{x}$ plane. At each iteration, I detect a number of peaks on the surface, each with an $x/\dot{x}$ measurement pair.

I need to select the peak which is closest to the current state as input for the next filter iteration. At the moment I am simply finding the difference in the $x$ and $\dot{x}$ dimensions, scaling them so they are similar in scale, and then computing the 2-D distance. This has been working fine for me, but I know it is not the "right" way to calculate the distance.

Is there a way to use the Kalman variables to determine which measurement is closest to the current estimate?

Not sure where this belongs, so I have this posted in stack, dsp, math, and stats exchanges.


texts: Transfer function of second order notch filter

texts:
Can anyone mention the transfer function of second order notch filter to remove the line frequency of 50 Hz, in terms of frequency and sampling rate.

Just like for Low pass Butterworth filter as,

$$
H= \frac{1}{\sqrt{1+\left(\frac{\omega_n}{\omega_c}\right)^4}},
$$

where $\omega_n$ is the signal frequency and $\omega_c$ the cutoff frequency.


texts: Convolutional coding that payload length does not change

texts:
This is the question about Sigfox radio specification.

The stack:

Specifically, in Section 3.10 of the specification, we have a convolutional coding (after CRC) with three possible polynominals.

I am wonder whether really error correcting code (ECC) because the output and the input lengths are identical (Figure 3-6 in the specification, or see below).

The question is that if this is ECC, could anyone tell me what is the code rate? And if this is not ECC, what is its role?


texts: Finding the phase and frequency of a tone with Doppler. What is wrong in the code?

texts:
I am trying to find the phase and frequency of a single tone having a small Doppler (10 Hz). I have the following code. The integration time for the operation (the length of the time vector) is defined in the code. It works if the integration time is a correct multiple otherwise not. Why is this failing?

Here is the code:

clearvars; close all; clc;

 Fc = 10.005e3;
 Fd = 10; % random doppler
 Fc2 = Fc+Fd;
 phi = deg2rad(37.8);
 disp(['True Freq: ',num2str(Fc2),' Phi: ',num2str(rad2deg(phi))]);
 Fs = 1e6;      % Max. Sampling Freq allowed
 Ts = 1/Fs;     % Sampling interval
 T_intg = 1.2347e-3; % Max. Integration time allowed

 t  = 0:Ts:T_intg;                        % Time Vector
 s  = exp(1j*(t*2*pi*Fc2+phi));           % Signal Vector

 Fs = 1/Ts;                               % Sampling Frequency
 Fn = Fs/2;                               % Nyquist Frequency
 L  = length(s);

 fts = fft(s, L)/L;                       % Normalised Fourier Transform
 Fv = linspace(0, 1, fix(L/2)+1)*Fn;      % Frequency Vector
 Iv = 1:length(Fv);                       % Index Vector

 amp_fts = abs(fts(Iv))*2;                % Spectrum Amplitude
 phs_fts = angle(fts(Iv));                % Spectrum Phase

 [val, index] = max(amp_fts);             % locating the peak in the spectrum
 rad2deg(angle(fts(index)))                % this gives 104 instead of 37 degrees


texts: Confusions regarding meaning of "kernel" in signal processing and image processing?

texts:
While studying filters in DSP and DIP, the term "kernel" is often encountered
What is the meaning of "kernel" in context of filters??Please kindly explain in simple words with an example

I have also attached a snapshot from gonzalez

What is meant by "neighborhood of operation" here?

and in last line,author says that mask,template & window terms are synonyms for kernel.

Is it true even in case of FIR windows such as Hanning window,Blackman window etc? Those windows are also kernels?


texts: Change duration/tempo of a .Wav file without introducing pops or pitch shift?

texts:
I am attempting to develop a similar function as Audacity's Effects>Change Tempo effect.

Or more specifically: I am in need of a way to create a function in my own project that can modify the length of a .wav file without changing it's pitch or introducing pops/leakage.

I am writing a standalone application (in C#) that already has some functions to create and save .wav files. But I just realized that I must add some processing to allow me to adjust the duration/length of a .wav so I can sync it to be exactly as long as an associated animation (which can not be altered).

I am very interested in knowing if there are any examples that exist which demonstrate concepts on how to change a .wav file's tempo (without noticeable distortion)?

I should specify that I would not need to adjust the length by much more than ±5%.

What I have tried already is to remove slices of data at regular intervals for a given "slice width", which kind of worked but introduced many pops as the sine waves were no longer smooth/continuous.

I have also tried simply adjusting the playback speed at the player by a multiplier to adjust to a desired duration, but this obviously cause the pitch to change (and somewhat noticeably too).

This is something that I am very anxious to understand as my project now seems to be in trouble unless I can somehow synchronize the audio. Could someone please help me figure out if I can write my own function to process this?


texts: How to measure correlation between 1 bit samples and theoretical original waveform?

texts:
I'm investigating quantization error.
I have an analogue waveform that looks like this and is of theoretically infinite resolution:-

I've sampled it as (8 bit oscilloscope readings & 0b1) to produce a file of single ones or zeros. The samples are partially random as per quantization theory, but I want to check the correlation with the analogue original. Theory suggests that there should be some correlation.
How can I do that? Clearly I can't compare the 1 bit samples with the 8 bit samples they came from.
Put another way, how can I show that my quantization is entirely uncorrelated to the underlying waveform?

Edit & maybe a solution:-
St is the theoretical and perfect(!) signal.
S1 is my sample set as (scope reading) & 0b1.
S2 is another independent sample set as (scope reading) & 0b1.
It occurs to me that I can't possibly compare S1 with St, as quantisation error will be present in both sets. But, if S1 ~ St, then by inference S2 ~ St. Hence I can test the hypothesis that S1 ~ S2 since they're both derived from common signal St. If there is a correlation, then S1 must be correlated to St too.
So I can just use normal correlation analysis (or Hamming distance, Jaccard index) between S1 and S2 to measure any. Yes/no?


texts: Design and a notch filter using two Butterworth filters in parallel (MATLAB)

texts:
I am doing a project to filter out a 1kHz noise inserted to a 30sec .wav audio file.
To filter out the 1kHz noise, I am required to design using two Butterworth filters connected in parallel where one is a low pass and the other is a high pass Butterworth Filter.

I have tried to notch out the 1kHz noise by setting the bandwidth to (+/-10Hz) of the 1kHz noise.
Therefore, the range of the noise signal to be attenuated will be from 990Hz to 1010Hz.

I am planning to obtain the transfer function in s-domain first and then transform it to z-domain.

The designed transfer function will be inputted into MATLAB and the filtered audio file will be obtained.

May I know the steps in designing this filter? Cost-effectiveness and ease of implementation of the filter is emphasised.

Thanks.


texts: Calculating an in-loop signal as part of a hierarchical control loop

texts:
I've got a control system with two feedback paths, with each path going to a different actuator that corrects the error in the system. One feedback path provides feedback at low frequencies, and the other at high frequencies. Here is the loop:

I want to calculate the signal on one of the actuators (the node at the output of the "slow" block) while the loop is closed, given a disturbance that I inject before the summing junction. The signal I want to calculate is "A". The input is normally zero (I just want to suppress disturbances to zero, not some set-point).

I initially started by calculating the open loop gain provided by the inner plant-slow-sum loop, i.e. A = slow * plant * disturbance, and then divided the disturbance by this open loop gain to give me something that approximates the disturbance suppressed by the closed loop, but this is not correct because in reality some of the disturbance is also suppressed by the outer plant-fast-sum-sum loop. In effect, the disturbance that the "slow" loop will see at high frequencies will be much lower, since the "fast" loop suppresses it there.

How do I calculate the closed-loop signal at A due to the disturbance, including the effect of the outer loop?


texts:
I'm reading about when a $H(z)$ system is min phase or max phase, and I get when a system is a min phase or max phase depending on the location of the zeroes. But I don't get at all the phase shift part. on one point wikipedia says that these two transfer functions:
$\frac{s+10}{s+5}  \; and \;\; \frac{s-10}{s+5} $
$\\\\$
"have equivalent magnitude responses but however the second system has a much larger contribution to the phase shift." So what is phase shift?


texts: How can I improve my motion capture system?

texts:
Please excuse me if this is not the proper place to ask this. I just need some general advice.

I have set up a very basic homemade motion capture system where I put three TrackIR cameras from NaturalPoint in a three fixed positions. I want to capture movement of a single object in three axis. So I laid out two of the cameras on both sides of my working area, plus the third camera above the object, pointing downwards.

I then use reflective markers to track the position of my object in the cameras.

I don't do any calculations on my system aside from a simple coordinate transformation from that of the real world to that of my virtual world. I do this separately, for each of the cameras. Since they are in fixed positions I already know how movements translate from one "world" to the other.

Is there any way that I can improve my mo-cap system? Should I apply some kind of triangulation (isn't triangulation implicit in my system already, given that Im using fixed positions for my cameras)? What filtering could I use to account for punctual lost of tracking elements?

Note that I want (ie. need) to stick to my current hardware.


texts:
I'm looking at the GNU Radio GMSK modulator block. The Gaussian taps are somehow convolved with a rectangular window with 4 * samples-per-symbol taps owing to the fact for GMSK, the filter spans four symbols. I was wondering if anyone could elaborate why do we convolve the Gaussian taps with the square wave?

self.gaussian_taps = filter.firdes.gaussian(
        1,             # gain
        samples_per_symbol,    # symbol_rate
        bt,            # bandwidth * symbol time
        ntaps                  # number of taps
        )
self.sqwave = (1,) * samples_per_symbol       # rectangular window
self.taps = numpy.convolve(numpy.array(self.gaussian_taps),numpy.array(self.sqwave))
self.gaussian_filter = filter.interp_fir_filter_fff(samples_per_symbol, self.taps)


texts: What is the advantage of Wiener filter for noise reduction of a 1D signal?

texts:
I have a 1D signal, acquired by an accelerator sensor that measure the vibration of a structure. What is the advantage of using a Wiener filter for noise reduction compared to the other (more classic) filtering methods?
I have a signal + noise $u[n]= s[n] +w[n]$ measurement and a noise only $w[n]$ measurement, taken at rest.


texts: Why is DTFT of $e^{jn\omega_0}$ an impulse train?

texts:
update : After asking the question, I figured out that DTFT result is an impulse train. Now my question evolved to, how it is derived in this way?

Using the DTFT formula seems not to be working, because of "sum does not converge".

$$e^{j\Omega_{0}n}=2\pi \sum_{k=-\infty}^{\infty} \delta\left(\Omega-\Omega_0-2\pi k\right)$$

However, that DTFT is defines as above in most places.

Why cannot I find the transform using the forward DTFT formula?


texts:
The coefficcients that define a FIR filter are $b_k = \{0, 0, 2, 0, −1, 0, 2\}$. The filter gets the following signal as input:
$$x[n] = −7\delta[n − 3]$$
Calculate the output of the filter (using the conv command), and draw the image of the input and output using the stem and subplot commands.


texts: The notation to express a datum filtered by arbitray number of one-dimensional filters?

texts: I'm a pure novice in signal processing but I'm writing a paper that needs to be expressed in the SP way. The question is, I want to express an N-dimensional datum goes through N one-dimensional filters, each along a different dimension of the datum. The simplest example is applying Gaussian blur filter on a 2-D image. You can separate the 2-D Gaussian filter into 2 1-D Gaussian filter and apply it along each dimension of the 2-D image one after the other. But my question is, is there an elegant notation to express such operation?


texts:
I need to build nonlinear model having recorded input/output data. The object is not simple and linear model will fail for sure.  What other options exist besides Least Mean Squares?

I need to identify dynamic system with one input. I’m provided recorded sequences of input and output data. I built linear regression model, which is called finite impulse response, but it is very inaccurate. I need more accurate model.


texts: Will lossy compression prevent me from performing audio analysis?

texts:
I'm starting a project where we will have to analyze bird recordings. We will try to measure the birds activity using the audio.

The recordings are really long (around 24 hours each) and recorded in IMA ADPCM WAV.
All files are about 1.5GB each.

When trying to import them into MATLAB, it simply locks up because of the size of the files (I think). So we were thinking about heavily compressing the files to 8KBPS mp3's, which would result in files of around 100MB.

Now I know mp3 compression is based on the psychoacoustics of the human hearing, so I was wondering whether using mp3 compression would throw away data dat would be needed for the analyzation.

If it does throw away data that we would need, what would be the best route to take? Are there other compression techniques that might be better suited to our needs?

Thanks in advance.


texts: Where did we get the DC term of the Accumulator from DTFT?

texts:
Define $y[n]:=\displaystyle\sum_{m=-\infty}^{n}x[m]$. The DTFT is found as follows:
\begin{align*}
y[n]&=\sum_{m=-\infty}^{n}x[m]\\ \\
&=\sum_{m=-\infty}^{n-1}x[m]+x[n]\\ \\
&=y[n-1]+x[n]\\ \\
\implies x[n]&=y[n]-y[n-1]
\end{align*}
Therefore, applying DTFT :
\begin{align*}
\mathcal{F}\{x\}[n]=X(e^{j\omega})&=\mathcal{F}\{y[n]-y[n-1]\}\\ \\
&=(1-e^{-j\omega})Y(e^{j\omega})
\end{align*}
Hence :
$$
Y(e^{j\omega})=\frac{1}{1-e^{-j\omega}}X(e^{j\omega})
$$
However, the correct answer is :
$$
\frac{1}{1-e^{-j\omega}}X(e^{j\omega})+\color{red}{\pi X(e^{j0})\sum_{k=-\infty}^{\infty}\delta(\omega-2k\pi)}
$$
My question is where did the DC term come from and why didn't it work in my derivation?


texts: Understading sampling theorem and aliasing

texts:
A real-valued analog signal with a flat spectrum between f = 0 Hz, and fmax is sampled at a sampling frequency fs = 24 kHz. The sampled signal is then processed with an ideal lowpass filter with cutoff frequency = 0.25pi.

How can we specify the range for omega in which aliasing can be torelated and also determine the maximum signal frequency fmax of the analog signal such that no aliasing components are present in the filter output?


texts: How would you generate this type of wave form?

texts:
I'm trying to generate a wave form that looks like this:

It looks to me like a square wave that has a large transitional section between going from "low" to "high" and back.

Any idea how to generate such a wave mathematically (and/or via some Python library)?

Note also that there is a higher-frequency carrier wave.


texts: Convolution and windowing using a buffer - how do I do overlap add?

texts:
I'm reading chunks of audio signal (1024 samples) using a buffer, applying a Hanning window, doing an FFT on it, then reading an Impulse Response, doing its FFT and then multiplying the two (convolution) and then taking that resulting signal back to the time domain using an IFFT. I'm zero padding the Impulse Response to match the length of the input signal chunk. However, I still get a noisy click sound after it reads every chunk, because of differences in amplitude I suspect.

At this stage I want to do an Overlap Add but am unclear on how to do it on a buffer. Convolving the input signal chunk with the Impulse Response results on a signal with the same number of samples (1024) so what am I overlapping when it gets to the next chunk?

Should I be using a longer impulse response and zero pad the signal
instead, and then add the rest of that signal after sample 1024 to
the next window and so on?  Or should I use a longer window on the
input signal?

I am doing all this in Python using PyAudio.

Thank you very much for your help


texts:
I know what a delay is, but
what's a delay line? What's the reason for calling it a "delay line"?


texts: determining FEC or parity bits from percent error

texts:
I am very new to signal processing and my background is in physics.  I would like to know if it is possible to determine the number of parity bits needs to theoretically get 100% transmission from a transmission channel with known percent error.  Over a given transmission 25-30% of the transmission is incorrect.

Is it possible to know with say a reed-Solomon technique how many paity bits I would need to add to each block to move in the neighborhoos of 100%?

I am processing my signal using MATLAB right now...


texts: Gaussian FIR filter with no multipliers?

texts:
What is the simplest way to implement a Gaussian FIR filter with unity gain coefficients and no multipliers?

Please preface your answer with spoiler notation by typing the following two characters first ">!"

Note: A Gaussian FIR filter is an FIR filter with an impulse response that is a Gaussian function.  By "unity-gain coefficients" I mean all coefficients in the filter structure are 1.


texts: Visualising OFDM subcarriers in the frequency domain

texts:
Im struggling to understand why in some images the frequency domain shows clear bars with amplitude and phase and with 0 frequency variation (it is clear that the frequency of the waves does not change in the image on the left).

Yet in others it shows it sloping or curving upwards:

Why isn't this a bar when the frequency of the signal clearly doesn't vary?

Same thing in OFDM. The constellation points should map to a single 312.5KHz aligned frequency and modify their amplitude and phase, yet when you look at the frequency domain on the OFDM carriers, they're sinc shaped waves instead of being bars. How is this possible?


texts: what is the best technique to detect a given signal : energy detection, autocorrelation, and matched filter based sensing, CRC?

texts:
In order to detect the presence of a signal impaired by noise, what is the best technique (best metrics) one can use ?

energy detection,
autocorrelation,
matched filter based sensing,
CRC
data aided or not ..
an explanation of the underlying algorithms would also be appreciated


texts:
Any help to get me started is appreciated. I have an array of values for the 4-QAM where attach is a copy of part of the values of the array. I want to know how do I get the bits 0 and 1s from those values. I would be using python to program. Thank you very much

I have searched the forum and the post below is something that I am trying to do..I understand on how to get the bits. Is there a formula where I can feed the above values?  Most references state to check the BER (bit error rate) but do we get the bits..

Recover IQ samples from signal

"If you just want to demodulate, then it is not necessary to recover I and Q components. What you have is most likely the QAM symbols. For example, for 16-QAM they are 16 different integer values that are assigned to different constellation points. You need to know the type of assignment (i.e. coding) that is used in the QAM constellation at the modulator."


texts: Suppose a signal is growing exponentially. We take its Z transform & find its ROC. So what are we supposed to do practically with the signals in ROC?

texts: Suppose a signal is growing exponentially. We take its Z transform & find its ROC. So what are we supposed to do practically? Are we supposed to manipulate our O/P by multiplying any signal from ROC such that our O/P converges & system is stable? And what if our signal is exponentially decaying so whats the purpose of finding ROC in this case (except for determing stability/casuality)?


texts: Confussion on roll-off factor, sampling time and Raised Cosine pulses.

texts:
I am reviewing some bandlimited communications theory and something is bothering me -

My book says increasing the roll-off factor makes the received sequence more tolerant to timing recovery errors but it's very clear that the strongest peaks are not at the sampling instances. Thus it doesn't seem that the system is very tolerant to timing errors regardless of the roll-off factor.

It looks like to make the peaks equate to the sampling times you would need to half the transmission rate.

Am I wrong? What's going on?


texts: Viterbi equalization of Fibonacci coding

texts:
I'm working a Viterbi equalizer for a BPSK signal which uses Fibonacci coding. In the more theoretical descriptions of the Viterbi algorithm I've read, there's a term for the probability of a particular observation from a given state. In the more practical descriptions I've seen this term does not exist. I assume it's because the probability of observing a 0 or a 1 are considered equally probable.

This isn't true with Fibonacci coding. I've read (but not seen the proof) that the probability of a 1-bit is ${1\over 2} - (1-\sqrt{1/5}) \approx 0.2764$.

Can this be used to advantage if the objective is minimizing the word error rate (words being delimited by two consecutive 1-bits)? Or would this make things worse, perhaps by making the detection of the word boundaries less likely to be detected, thus increasing the probability of word errors?

And if so, within the context of channel equalization or decoding, the Viterbi algorithm is used to find the minimum path metric through the trellis, with the path metric being the hamming or euclidean distance. More generally though, the Viterbi algorithm is used to find the maximum path through the trellis, using the probability of each transition. How to reconcile these two contradictory concepts?


texts:
I'm reading through the literature for finding the best way to calculate SNR for a ECG signal. While there are many algorithms for improving SNR, I couldn't find any that details any standard procedure for evaluating the SNR itself. The following are my questions.

What is the SNR of a typical ECG signal from the MIT/BIH database? What is considered a good SNR for a typical ECG signal?
http://www.physionet.org/physiobank/database/mitdb/
What would qualify as Signal and Noise for ECG? The assumption that I worked with is to apply a bandpass filter with cut-offs fc=[0.5Hz,45Hz] on the recorded ECG signal to get the Clean ECG signal, then subtract the clean ECG signal from recorded signal to obtain noise and then use the matlab snr function to evaluate  SNR with the clean signal and noise as arguments. But, in some papers it is recommended to reject everything except the QRS complex to get the clean signal. Is there any standard/recommended procedure for calculating SNR for ECG?
http://en.wikipedia.org/wiki/QRS_complex
What would be the best way to calculate ECG SNR in real-time in an embedded platform?